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1 Introduction 


High order finite difference methods (HOFDM) are ideally suited for simulations of complex physics requiring high 
fidelity solutions, and for a simple reason: efficiency. Low-order techniques (whether they be finite difference, finite 
volume, or finite element) lack sufficient accuracy for a given cost to make them competitive in this environment. 
HOFDM have for decades been used for direct numerical and large eddy simulations (DNS/LES) of fundamental 
problems in fluid mechanics. In recent years, they have begun to impact a wide variety of other aerodynamic ap- 
plications, including 1) subsonic noise reduction, 2) air space utilization/safety, 3) uninhabited combat air vehicles 
(UCAV), 4) circulation control, 5) fluid/structure interactions, 6) turbomachinary, 7) aircraft in flight, 8) helicopters, 
9) atmospheric chemistry, 10) windmill optimization. Other disciplines that HOFDM are making inroads into include 
11) computational electromagnetics, 12) ocean models, 13) meteorology, 14) magnetohydrodynamics, 15) weather fore- 
casting, 16) oil recovery simulation, 17) modeling of shallow waters, 18) transport of contaminant in porous media, 
19) viscoelastic flows, 20) semiconductor device simulation, to name a few. 

The basis of the argument used to advocate the use of HOFDM is not new, and in fact dates back to early 1970’s 
work by Kriess and Oliger [25] . They showed that the efficiency of HOFDM increased with each consecutive order up 
to approximately 5th-order for time-dependent hyperbolic problems. In spite of their theoretical advantage, however, 
two significant impediments continue to face HOFDM. They are: 1) stable boundary closure formulae maintaining the 
formal accuracy of the interior scheme, and 2) extension to complex geometries. 

Many theoretical (and ad hoc) boundary closures for HOFDM have been proposed over the past three decades. 
Most, however, still rely on low-order boundary closures to achieve stability and robustness. Unfortunately, even one 
low-order point can have a detrimental effect on asymptotic solution accuracy. Specifically, Gustafsson [18] showed 
that for hyperbolic equations, boundary stencil accuracy can differ from interior stencil accuracy by one order at most, 
if design order is to be maintained (e.g. (p) t/l -order interior with (p — l)* ,! -order boundaries maintains (p) tft- -order 
asymptotic accuracy). Svard and Nordstrom [38] recently showed how to relax this constraint to allow (p— 2) th - 
order boundary closures for strictly parabolic problems. Nevertheless, numerical schemes for hyperbolic, incompletely 
parabolic (Navier Stokes: NS) or parabolic equations all suffer the potential consequence of low-order closures; and 
thus, boundary closures must be carefully designed. 

High-order biased stencils near boundaries (on uniform grids) are susceptible to Runge-type oscillations/instabilities 
that frequently compromise the stability of the simulation. A simple remedy to the problem of boundary closure stencils 
originates in the early work of Kriess and Scherer [26]. They realized that stability issues could be controlled by simply 
constructing semi-discrete finite difference operators that satisfy a semi-discrete energy estimate. That is, the semi- 
discrete operator ” mimics” the discrete operator in terms of the energy of the system. The critical element that enables 
the existence of the continuous energy estimate is the integration-by-parts property of the differential operator. Thus, 
semi-discrete operators that mimic this property are said to satisfy a summation-by-parts (SBP) condition. Although 
Kriess and Scherer focused exclusively on low-order methods, Strand [36, 37] in his Ph.D. thesis and Olsson [32, 33] 
extended these ideas to HOFDM, and the field of summation- by-parts (SBP) high-order schemes was born. 

In principle, SBP operators can be derived for virtually every conceivable operator type. By construction, weak 
form FEM operators automatically satisfy the SBP property. Additionally, operators that satisfy the SBP property 
exist for 2 4 th -, 6 th - and 8 t,l -order central schemes [36, 10], Pade schemes [10, 1], upwind schemes [27], dispersion 
preserving schemes [29], and strong form nodal ’’spectral element” schemes [11, 15, 12, 20, 21, 22, 23]. Operators that 
satisfy the SBP form need not only be applied to hyperbolic systems. Indeed, parabolic equations are amenable to 
SBP forms, either by the action of two first order SBP operators or by constructing operators specific to the second 
derivative [10, 27]. Special considerations must be accounted for when dealing with these equations, but no great 
difficulties are encountered [28]. 

An implementation detail related to SBP operators emerged in the work of Carpenter et al. [10], dealing with 
the imposition of the physical boundary data. Specifically, they noted that if the physical data were ’’injected” in 
the conventional manner, then the energy of the semi-discrete system could become unbounded. Borrowing ideas 
from the interior penalty approaches used in Finite-Elements (IPFEM) [16, 4, 39] and those used in pseudospectral 
collocation [17], they devised a means of imposing the data, [which was designated SAT Simultaneous Approximation 
Term (SAT)], that guaranteed bounded energy while maintaining formal accuracy. The SAT procedure solves a linear 
combination of the boundary conditions and the differential equations near the boundary. A similar alternative was 
developed by Olsson [32, 33] based on orthogonal projection techniques, that also guaranteed stability and accuracy 
of the discretized system, including physical boundary data. 
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The extension of HOFDM to complex domains proves more difficult to accomplish. Multi-domain HOFDM are 
a natural candidate for extension, but are fraught with the complexity of interface coupling procedures. Specifically, 
they must address the interface accuracy constraints, as well as contend with coupled stability conditions and interface 
conservation issues. Many ad hoc procedures have been developed to accommodate multi-domain methods, including 
conforming, nonconforming, and overset/overlapping procedures. Once again, the key to success for these procedures 
seems to be reduction of the order of accuracy at the interface region and the subsequent loss of solution accuracy 
throughout the entire domain. 

A simple solution to the multi-domain problem for HOFDM was provided in the work of Carpenter, Nordstrom, 
Gottlieb[13] (henceforth referred to as CNG). Motivated by the recent popularity/success of Discontinuous Galerkin 
Finite Element methods (DGFEM), and the internal penalty approaches, they proposed combining the SBP operators 
within each domain, with a penalty technique at the interfaces. The formulation requires only weak grid continuity at 
the interface (Go: matching but not necessarily smooth), and is applicable to any operator that satisfies a SBP property. 
Furthermore, the operators need not be the same in adjoining domains. It maintains formal accuracy through the 
interface, is conservative and maintains a bounded semi-discrete energy estimate consistent with the underlying SBP 
operators. Although originally developed for the linear, constant coefficient, scalar, advection-diffusion operator on 
multiple domains, the CNG interface formulation was subsequently generalized to systems of equations with constant 
coefficients by Nordstrom et al. [30] , and finally to the linearized Navier-Stokes equations. [31] Similar developments 
within the nodal spectral element community were reported by Hesthaven et al. [22, 23, 24], 

In summary, the CNG formulation, patterned after the IPFEM/DGFEM methods, attempted to provide a ’’general” 
technique to couple multiple SBP domains. Conventional penalties were constructed at interfaces based on Dirichlet 
and Neumann data from either side. After accuracy, stability and conservation were enforced, two free parameters 
remained, one for convection and one for diffusion. Significant maturation has occurred within the DGFEM community 
over the past 6 years, especially in the treatment of diffusion terms. Recent work by Arnold et al. [3] provides a clean 
unifying description of all the approaches typically used in IPFEM and DGFEM, and has tied up many of the loose 
ends. These, and other recent developments for the treatment of the DG diffusive terms, make it apparent that the 
CNG formulation is (although correct) not a general procedure. In fact, additional penalty terms could have been 
included in the original formulation! 

The objective of this paper, is to present a ’’general” treatment of multi-block interfaces, connecting adjoining 
domains discretized with SBP operators. The strong (nodal) form of the governing equations is the primary focus of 
this work, although an elementary discussion of FEM methods is presented to motivate the new interface treatments. 
Two equivalent formulations are presented for the new interface conditions, similar in spirit to the flux and primal 
forms used in the FEM literature. Both new formulations are valid for the linear, constant coefficient, scalar advection 
diffusion system, and are accurate, stable and conservative, just as was the previous CNG [13] formulation. Five free 
parameters remain unspecified at interfaces in the new formulations, and can be chosen quite arbitrarily (only mild 
stability constraints). Furthermore, judicious choices for the parameters yield commonly used schemes from the FEM 
literature. 

The rest of this paper proceeds as follows. In section 2, we present the unified DGFEM approach to the interface 
treatment, written in the weak form. This formulation is then reformulated in terms of the strong form of the 
equations. In section 3, the new approach is presented for the strong form of the equations, including proofs of 
stability and conservation. This approach incorporates penalties on the solution and first derivative that are more 
general than the original CNG formulation. In section 4, a general penalty procedure is introduced that includes 
penalties on the second derivatives. Section 5 presents a derivation of the theoretical accuracy of the new penalties, in 
the context of HOFDM. Section 6 identifies several commonly used schemes in terms of the new penalty parameters. 
In section 7, numerical examples are presented for a variety of SBP operators. A discussion followed by conclusions 
are presented in sections 8 and 9, respectively. Finally, several appendices include related theoretical derivations. 


2 Discontinuous Galerkin FEM 

Our goal is to develop accurate, stable and conservative interface treatments applicable to multi-block HOFDM and 
other general collocation approaches; methods formulated in terms of derivatives of the strong form of the equations. 
Given the close similarity between strong form HOFDM and weak form DGFEM approaches, an elementary discussion 
of discontinuous finite elements for parabolic equations is presented. The derivation will in turn motivate a new penalty 
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treatment for the strong form of the equations. 


2.1 Nomenclature and definitions 

Parallel development of two distinct methodologies and nomenclatures exists in the discontinuous finite element litera- 
ture, devoted to solving elliptic and parabolic equations. The biharmonic equations are the genesis for Interior Penalty 
(IP) methods, which date back as early as the middle 70’s. The IP methods (and others based on the biharmonic 
equations) are often described as using the primal formulation , and they are primarily used for solving elliptic equa- 
tions. Interest in parabolic equations with a strong hyperbolic component has recently motivated the development 
of flux formulation methods, named for their utilization/reliance on inter-element flux reconstruction techniques such 
as exact or approximate Reimann solvers. Arnold et al. [2, 3] present an extensive review that unifies nearly all the 
DG schemes that appear in the literature. They study both flux and primal formulations, establishing the precise 
relationship between the two approaches. 

A brief derivation of the primal formulation and flux formulation is now presented. The objective is to establish 
the relationship between these two approaches, and a third methodology; the strong formulation, which will be the 
exclusive focus this work. We begin by discussing the flux formulation, and then, derive the primal form and strong 
form as a special cases of the flux formulation. 

Consider the scalar, 1-D, parabolic equation, 


u t = u xx , \x\ < l,t > 0. (1) 

with suitable initial data and boundary conditions. (For ease of presentation, but without loss of generality, we study 
the 1-D, linear, strictly parabolic equation. ) An example of a flux formulation FEM is the one proposed by Cockburn 
and Shu [14], denoted the local discontinuous Galerkin (LDG) method. Their approach (modeled after that of Bassi 
and Rebay [5]) begins by transforming the second-order equation to a ’’system” of first-order equations 

u t = dx, q = u x , \x\ < l,t > 0. (2) 

Multiplying equation (2) by test functions, and integrating by parts, the LDG method is then defined as: “find 
u, q € Va x such that for all test functions v, w € Vax> the following relations hold” 

f IjUt vdx = ~ fjj qv x dx + q j+ ±v7 +i - qj-xv+x 

f T qwdx = — fr uw x dx + u i4 _iw~, i — ft,- iu> + i . w 

with q and u being suitably chosen interface fluxes. Figure 1 provides a schematic of the nomenclature for the 

discretization. Here, the intervals Ij = [xj_x,Xj + x] is defined for j = 1, ...,1V, on the periodic domain H = [0, 27 t] . 

The points Xj are defined by Xj = \{Xj_ i + Xj + x ) and the element spacing by Ax.j = Xj + x — Xj_ i . The test 
functions v,w are piecewise polynomials of at most degree k such that v, w € VAx, where 

VAx = : v is a polynomial of degree at most k for x £ Ij,j = 1 (4) 

Define ej to be the edge between elements j and j + 1, and the union of all edges to be ej = T. Further define 
the vectors nlj and n2j on ej, pointing exterior to each element adjoining ej, and the superscripts “ and + to be the 
limiting interface value approaching from the left and right, respectively. 

Next, note that both the solutions u, q and the test functions v, w are discontinuous at the element interfaces. To 
facilitate manipulation of equation (3), we define average and jump operators on the edges. Written in terms of an 
arbitrary variable these operators are defined as 

K\ ~ k ' -V IK*]] = (C^+CM) Ve ; ( 5 ) 

where (ff are the left and right values of the double valued function Q on the interface j. Summing equation (3) over 
all elements defines the method on the entire domain Q. 

f Q u t vdx = - f Q q v x dx + X)r Hv w ]] 

f n qwdx = — / 0 uw x dx + JDr [Ml 

4 


( 6 ) 
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Figure 1: Two elements showing DGFEM nomenclature. 


e j + 1 


The primal formulation can be obtained from the flux formulation as follows. Allowing the test function w in 
equation (3.b) to satisfy w = v x (assuming v ^ cq, with cq a constant) produces the equation 

Jq q v x dx = — f Q uv xx dx + [[wfa;]] (7) 

which when substituted back into equation (3. a) yields 

fn u t . v dx = f Q uv xx dx - )Cr( [[«w]] - [[?«]]) (8) 

Integrating equation (8) by parts once more, produces the primal formulation 

fn ut. v dx = — f Q u x v x dx - X)r( [[(* - «)«®]] - [[&]]) ■ (9) 

The exact form of interface fluxes u or q is immaterial, for the purpose of the previous derivation of the flux and 
primal formulations. In practice, however, these fluxes can not be chosen arbitrarily. Additional constraints (such as 
global stability and conservation), establish admissible relationships between the u or q. See Appendix A for a brief 
discussion on admissible fluxes satisfying global stability, as well as a listing of several popular schemes. 


2.2 Strong Form 

The previously discussed methods, are defined based on integrals of the weak form of the governing equations. Specif- 
ically, equation (1) is multiplied by an arbitrary test function, and then formally integrated once or twice. We now 
derive a strong form equation suitable for developing collocation approaches. 

Starting with equation (9), an auxiliary primal form can be obtained by once more integrating the term u x v x . 
resulting in the expression 

fn u t v dx = f n u xx v - [[(^ - u)v x ]\ - [[(q - u^]]) • ( 10 ) 


Reverse engineering strong form methods is most easily accomplished at the element level. Reformulating equation 
(10) for a single element yields 


fj Ut v dx = Jj u xx v dx 


[(“-“);+!%•+! - % + _i] 

i(Q ~ u x )J +i vT + t - (q-u x )+_ 


( 11 ) 


which upon substitution of the primal flux form given in equation (153) (and reversing the order of multiplication), 
yields 


fj v Ut dx = fj v u xx dx 


%.+ i(Cl [[u}\ + C 2 M) i+ 1 
wt i(ci [M] + c 2 [[«jti 

J 2 J 2 

ViMM] + c 4 [[u x ]])~ + 1 
^_i( c 3[M] + C 4 [[uJ])t i 

J 2 J 2 


(12) 


The inverse transformation from weak to strong form involves replacing all polynomial integrals and derivatives, with 
discrete matrix vector operations on collocated data. The only subtle aspect of this task is representing the interface 
data in collocation form. To accomplish this, we define the discrete boundary operators 


e+ = [0, 0, ..., 0, 1] 


= M,...,o,of 


5 


(13) 



The action of these vectors isolates data at the left (— ) and right (+) edges of the element. 
Next, assume that discrete derivative and integration operations are defined by 


( x = V( + 0(h r y, J^vCdx = r, T V C + 0(h r ) 


( 14 ) 


where T> and V are discrete nodal differentiation and integration operators, respectively. 

Expressing equation (12) in terms of these definitions (13), (14) and the element identity matrix T.j , yields 

Vj T { ?U f = POTU + Z> T e+( Cl [[U]] + C2[[VU\]) j+ i 


V T e-( Cl [[lJ}] + C2[[2>U]]). _i 


X T e+(c 3 [[U]] + c 4 [[D U]]). +i 
^eJ(c 3 [[\J}] 


>o+\ 
C4[[2>U]]),-_1 


(15) 


} 


Equation (15) is satisfied for an arbitrary test function V T . A sufficient condition to satisfy the equation is for the 
bracketed terms to be identically zero. Dropping the premultiplication by the test function V T results in the strong 
form representation of the Discontinuous Galerkin FEM. 

The HOFDM and the nodal spectral element communities are accustomed to formulations (and stability proofs) 
involving two adjoining elements and the included interface. Rewriting equation (15) in this form yields the desired 
result 


T„+ 


V 
l T e+ 


VV i +l 


= vvvu j+1 


J 

D T e 

rT 


(aim + c2[[d u]]) i+ i 
tteim + c 4 [pu]]). +l 

:- + l( Cl [[U]] + C 2 [[D U]]). + | 
:- +1 ( C3 [[U]] + c 4 [[DU]]) j+ . 


(16) 


Note that the derivative and integration operators T> and V remain undefined in equation (16), so that this method 
could be applied to a wide variety of difference operators. Clearly, the method defined in equation (16) mimics the 
properties of the DGFEM method defined in equation (10). It is, however, legitimate to question whether the two 
forms are equivalent. 

Consider an arbitrary set of collocation points on the interval [xq = x j-±i x h x 2, x n- i , x n = x j+ 1 ], and let 
f(x) be defined everywhere on the interval. The interpolation polynomial / at (: e ) that collocates f(x) at the points Xj 
is given by 


N 


f N (x) = JnJ = Y^ f(xj) L j(x) . 

j= 0 


(17) 


where Lj(x) are the Lagrange polynomials defined on the discrete set of points. It is well known (see Carpenter and 
Gottlieb [12]) that for this case, 

f* = Z>f (18) 

where the differentiation matrix T) is given by 


V = V 


-l 


e = ^L k ( Xj ) 


with 


ri+ 2 


Pki = 


Li(x)L k (x)dx ; q k i = 


'o-i 


p+2 Q_ 

dx 


Li(x) L k (x) dx 


(19) 


( 20 ) 


'o-i 


Substituting these definitions back into the strong form given by equation (16), and comparing the resulting expression 
with that of equation (10) (using the Lagrange polynomials as test functions), we achieve term by term agreement 
with the original weak form method. Thus, the strong and weak forms of the linear, constant coefficient heat equation, 
result in identical solutions. 

Remark. One has only to include a variable coefficient in the equation for the two formulations to produce different 
results. In this case, the collocation approach would have a truncation error associated with the inexactness of the 
differentiation operation. The DGFEM would, however, integrate exactly the polynomial of higher degree. 
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Closing Remark. 

The principle result of this brief review of DGFEM is equation (16), which expresses the “strong form equivalent” of 
the underlying FEM method. In this expression, it is clear that the V T matrix plays a critical role in the collocation 
approach, in determining accuracy, stability and conservation at the interface. The previous derivation of “general” 
multi-domain interface conditions for HOFDM, (see Carpenter et al. [13]), although correct, lacked the generality 
necessary to include the V T interface contribution! Armed with the V T matrix in our arsenal of mathematical 
tricks, we now derive the general interface coupling conditions for arbitrary collocation approaches, including HOFDM 
and nodal spectral collocation methods. 

3 The semi-discrete problem 

Our goal is to derive stable, conservative and accurate semi-discrete interface coupling operators. As we shall see, 
these operators closely relate to (mimic) the corresponding properties of the continuous operators. In the discussion 
to follow, we shall introduce each step with a discussion at the continuous level, followed by the derivation of the 
semi-discrete mimetic operators. 

3.1 Definitions 

Consider the linear initial boundary value problem 

Wt = Pw + SF(x,t ) ,x e!l , t > 0, 

w = Sf{x ) ,xGTl , i = 0, (21) 

Lqw = Sg(t) , rrST , f>0, 

where P is the differential operator and Lq is the boundary operator. The initial function 8f, the forcing function SF, 
and the boundary data Sg are the data of the problem; w denotes the difference between a solution with data f,F,g 
and one with data / + Sf,F + Sf,g + Sg. There are many concepts of wellposedness, see [19]. Here we consider the 
following definition. 

Definition 1 The problem (21) is strongly wellposed if the solution w exists, is unique, and satisfies 

IMIn + f IMIr dt < K c e^{\\Sf\\l + f{\\SF\\l + \\Sg\\ 2 r )dt}, (22) 

Jo Jo 

where K c and rj c may not depend on 5F, Sf,6g. || • \\q and || • ||r are suitable continuous norms , such as 

/•+ 1 p+i 

\\uv\\n = J uvdx, |M|n = y uudx, ||w||r = Ma=-i + M*=+i 

The senhdiscrete version of (21) is 

( Wj)t = Qwj+5Fj(t) ,XjGtl ,t> 0, 

Wj = Sfj ,ij e!l ,f = 0,, (23) 

LpWj = 5g(t) ,XjGT ,t> 0, 

where <5 is the difference operator approximating the differential operator P, SFj is the forcing function, Sfj the initial 
function, Lp the discrete boundary operator where numerical boundary conditions are included, and Sg the boundary 
data. It is assumed that (23) is a consistent approximation of (21). 

Closely related to the concept of wellposedness is the concept of stability. 

Definition 2 The problem (23) is strongly stable, if for a sufficiently fine mesh, the solution Wj satisfies 

\H\ 2 v + f IMIr dt < K d e^{\\5f\\l + f\\\6Ff r + \\Sg\\ 2 )dt}, (24) 

Jo Jo 

where Kd and rjd may not depend on SFj , 5fj ,Sg. || • \\j> and || • ||p are suitable discrete norms, such as 

\\vw\\ v = v T Vw, \\u\\ 2 v = u T Vu, IMIr = M*=-i + Mx=+i. 

Here, ||rw||-p, ||u||p, and ||tt||p denote the L 2 scalar product, the L 2 norm, and the boundary norm, respectively. 
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3.2 Discrete Spatial Operators 


Let U and T>JJ be the discrete approximations of the scalar quantities u and u x . The approximation T>U of the first 
derivative 

VU = QU, "Pu x - Qu = TT ell |T el | = 0(Ax i ,Ax b ) (25) 

satisfies the SBP rule 


(U,Z>V)„ = U N V N -U 0 V 0 -(VU,V) v (26) 


(U, V)-p = U T W, V = V T , Q+Q t = B , B = diag[- 1,0,. ,0,1], (27) 

0 < Pmm Ax||T|| < ||'P|| < p ma a;Ax||X||. The matrix V is symmetric positive definite, the matrix Q is skew-symmetric 
with the exception of boundary points. 

A second derivative operator can be obtained by applying any two first derivative operators, or by constructing 
the operator directly. Both techniques are suitable for spectral collocation. In HOFDM, using two consecutive non- 
dissipative first derivatives, results in a second derivative operator that is unnecessarily wide and inaccurate and can 
lead to odd-even mode decoupling. The operator, however, often leads to stability if varying coefficients are considered. 
Forming the HOFDM second derivative directly requires an operator with the following properties, 

V 2 U = p-^U, -Pu xx - 77 2 u = VT e2 , T e2 = 0{ Ax\ Ax b ), (28) 

77 2 = {-S t M + B)S, (29) 

as suggested in reference [13] (see also [28]). The matrix B is given in (27); M is positive definite, i.e. , U T A4U > 0 
and 0 < m m i n Ax\\l\\ < ||jw|| < m max Ax||X||. 

The matrix S is nearly diagonal with a discrete representation of the first derivative on the first and last rows, 
{Su}o = {Vu} 0 = u x (x 0 ,t) +T e 3 , {Su}n = {Vu}n = U x (x n ,t) + T e3 , 
where |r e 3 1 = 0( Ax r ) and 

soo Sol S 02 S03 

0 1 0 

0 1 0 

0 1 0 
0 1 0 

Snn—3 Snn—2 Snn — 1 Snn 

The second derivative defined in (28) and (29) satisfies a modified SBP rule We have 

(V,V 2 V) v = U„{DV}„ - Uo{PV}o - (SVfM(SV). 

The notation |T e i|,|T e 2 | = 0( Ax l ,Ax b ) and |T e 3 | = 0(Ax r ) means that the approximation of the differential 
operator is accurate to order i in the interior of the domain, to order b at the boundary and that the approximation of 
the boundary conditions is accurate to order r. This nomenclature is valid for any of the SBP schemes mentioned in 
the introduction, including HOFDM and spectral collocation. The relation between the different orders of accuracy, 
i.e., i, b , r is discussed in section (4) below, where formal proofs of accuracy are given for HOFDM. 
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3.3 The semi-discrete single domain problem 


The following semi-discrete derivations and proofs are closely related to proofs of wellposedness at the continuous level. 
For completeness, a complete derivation of wellposedness is presented in Appendix (B) for both a single and a two 
domain problem. A derivation of the semi-discrete equation is now presented. 

Imagine that we are solving the constant coefficient, linear Burgers’ equation 


Ct + a c x 

= cc xx + F(x,t) 

O 

Al 

-to 

,-l < X < 1, 

c(x, 0) 

= f(x) 

0 

II 

•So 

,-1 < X < 1, 

7 C (— 1, t) ~ ec x (-l,t) 

= 9-i{t) = i-i(c) 

O 

Al 

-to 

,X = -1, 

C c (+1, t) + cc x (+l, t) 

= g+ i(t) = L +1 (c) 

0 

Al 

-to 

,x = +1, 


with constants a and e (0 < e). The boundary conditions that lead to a wellposed problem satisfy the following 
constraints [(a detailed discussion of wellposedness is included in Appendix (A)]. 

0 < a + 2 £ ; 0 < -a + 2q . (31) 


Consider the semi-discrete, finite difference representation of the linear Burgers’ equation (30) of the form 

C t + a VC = eVVC + F(C) 

- V~\a - 1 [LB X (C) - g_ x ] e_!} 

- V~ 1 {a + 1 [L+! (C) - 5 +i]e+i} ,t > 0 -1 < x < +1 { ’ 

C(a+0) = i{xi) ,t = 0 — 1 < x < +1 . 

The discrete boundary operators are defined by 

i-r(C) = (7 C- c2?C)|_ 1 ; L^( C) = ((C + eZ>C)| +1 , (33) 

while the penalty vectors are defined by e_i = [1, 0, ..., 0, 0] T and e+i = [0, 0, ..., 0, 1] T . 

We have introduced a uniform mesh —1 < Xi < +1 with Xq = — 1, x n = +1 and 27 = —1 +iAx, and the boundary 
conditions are treated weakly via a penalty technique (see reference [10] ). 


Theorem 1 The approximation (32) of the problem (30) is strongly stable if the boundary conditions satisfy the 
constraints given in equation (31), and the numerical boundary conditions given in equation (32,33) satisfy the following 
inequalities: 


with 


1 _|_ 



< cr_! 


< 1 + 222 
— e 




a 


1 

ei_ T V ei 


(34) 


Proof : The energy method applied to equation (32) leads to 

j- t \\C\\ 2 v + aC T [VV + V T V]C = 2eC T [VVV + V T V T V]C + 2||C,F|| p 

- 2a- 1 C- 1 [LB 1 {C) - g-i] (35) 

- 2 a +1 C +1 [L^(C) - g+i]. 

The definitions of the first derivative operators V = V~ 1 Q, presented in equations (25) - (27) lead to 

C T (PV + V T V}C = C t BC 

c T [rvv + v T v T v]c = c t bvc - c T v T vvc 

Introducing equations (36) and (33) into (35) yields 

&\\C\$ = —2 e (VC) T VVC + 2||C, F||p 
+[— aC T £>C + 2 e C T VC}\ X X Z±{ 

—2 U-! C-i [(7 C- ePC)U - 5 _r] 

— 2 0+1 C + \ [(CC + eVC)\ +1 — 5+1] 
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(37) 



Substituting the relations 


2 II C, F||p 

2C'_ig_i 

2C+15+1 


-|IV5?C - + 77|| C HI, + i||F||2 < 

~[y/n ^ c - 1 - 9-if + v- (C'-i) 2 + ^ (g-i) 2 < 

^[\/V+C' + i - -^=<?+l] + ?y + (C + 1) 2 + (ff+l) 2 < 


+ ?? ||C||2 + I||F||2 

V^( c -i) 2 + 7b(^-i) 2 
V^( c + 1 ) 2 + tIt(5+ 1 ) 2 


into equation (37) yields the equation 

eVt jkie-^WCWv} + 2 e||PC||| 


+ 7,\\n 2 r 

\WvC-^Ff v 

TT (9- 1) 2 

- <T-1 [0TC_1 - ^=3-i] 2 

^ (5 + i) 2 

- cr+1 [y/d+C+l - T^yff+l]" 

$-1 

$+1 


with 


C- 1 

T 

—a + cr_i (27 — rj-i) +e(l - cr_i) 

C-l 

VC\-i _ 


+e(l — fj_ x ) +2 act 

_2>G|_ 1 _ 

' C +1 ■ 

T 

+a + ct+i( 2C — ?j+i) — e(l - 0 + 1 ) 

' c +1 ' 

vc\ +1 _ 


— e(l — cr+i) +2ea 

DC\ +1 _ 


( 38 ) 


(39) 


(40) 

(41) 


It is shown in Appendix (B) that part of the quadratic diffusion term 2 e ||2?C|jp can be moved to the (RHS) to 
stabilize the boundary terms $_i and 4?+i. Only a small portion can be borrowed, however, without destroying the 
definiteness of the term. The new diffusion term remains definite in the norm 


V = V — a e ; _ ei_ T 


See Appendix (C) for the value of alpha that produces a semi-definite matrix V . 

To ensure stability the terms $ + i and $_i must be positive definite, and thus lead to stability if the inequalities 
given by equation (34) in Theorem 1 are satisfied. (Here we assume that the parameters r/±i = 0.) Note that the 
wcllposedness conditions given in equation (31) ensure positivity of the square roots in both expressions. 

Grouping like terms results in 


Wt{e~ Vt \\C\\ 2 v} + e ~ vt {2 e \\VC\\ 2 p + BC d } = e~ * {Data - PosDef d } 

with 

BC d = $_i + $+i 

Data = +±\\F\\ 2 v + ^g +1 2 + ^g^ 2 

PosDef d = +|| - ^F\\ 2 v + a +1 (0£T C +1 - ^g+if + <t_i (0TTC - 1 - ^=9- if 


(42) 


(43) 


Finally, integrating equation (42) in time leads to 

\\C{x u T)\\ 2 v + / 0 T e-^- T ) {2e\\VC\\% + BC d }dt = e+* r {||/(* i )|£ + e " * [Data - PosDef d \dt} (44) 

from which follows a strongly stable estimate of the form 

\\C\ \ 2 v + [' \\C\\ 2 r dt < K d e^{\\f\\ 2 v + f{\\F\\v + IMIrM, (45) 

Jo Jo 

Hence (32) is a strongly stable approximation. □ 
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Xn — X 0 


+ 


Ax l i Ax r 


H 1 1 1 1 1 1 1 1 h 

Xn — 5 X n —3 X n — 1 1 X\ X2 X3 

I 

x = 0 


Figure 2: The mesh close to the interface at x = 0. 


3.4 The semi-discrete two domain problem 

3.4.1 A ’’Primal” discretization technique 


We are interested in the stability, accuracy and conservation of the new interface treatments. The previous single- 
domain result established the conditions required to achieve semi-discrete strong-stability, subject to an arbitrary 
source term and consistent, wellposed initial and boundary conditions. For clarity of presentation (and without loss of 
generality), we shall neglect the source term, and the physical boundary conditions from the subsequent two-donrain 
derivations. 

A collocation approximation of the coupled problem (168) written in primal form and motivated by the FEM 
proposed by Baumann and Oden [7] is 


Ti\5 t + aTiVW - eTiVW 

= 

loo[Ui~Vi\ 

Ii ei_ 

+ 

loi[(PlU)i - (‘DrV)i] 

Vi ei_ 



+ 

ho[Ui — Vi] 

Vfe,_ 

+ 

himu), - (D r V)i] 



U(s,0) 

l~*r Vj + a'Pr'DrV — cT 

= 

0 . 

roolVi-Ui] 

T, 

+ 

xoi[(V r V)i — (D r U)i\ 

T, ei + 

(46) 


+ 

nolVi-Ui] 

Vr ei + 

+ 

ru[(V r V)i - {V r U)i\ 

Vfe i+ 


V(Z,0) 

= 

0 . 






The variables in the left (subscript 1) \xq = 

-1, 

x n = 0 ] and right (subscript 

r) [x 0 =0,x m = + 1 ] 

domains are 

U and 


V, respectively, (see Figure 2). 


Definitions of the ni ( n r ) dimensional operator T>i ( T > r ) are given by 


Vi = V l 1 Qi\ -1 < x < Xi ; V r = V r l Q r ; Xi < x < 1 

and e;_ and ei + are given by e*_ = [0, 0, ..., 0, 1] T and e* + = [1, 0, ..., 0, 0] T . Note that the difference operators T>i, 
and T> r are in general different in the left and right domains (A xl 7 ^ Axr and/or iii 7 ^ n r ). 

What remains is to determine the values of the interface parameters loo, loi, ho, hi, and r 00 , roi, rio, rn, that 
lead to stability and conservation. 


3.4.2 Stability 

Remark. Stability of the one domain problem does not imply stability of the multiple domain problem. Stability 
means that the solution can be estimated in terms of the (bounded) boundary data. In a multiple domain problem, 
the boundary data are made up of the solution(s) in the other domain(s). Boundedness of the data would require an 
a priori assumption. 

Theorem 2 The approximation (46) of the problem (30) is strongly stable if the eight parameters are related by the 
two equalities 

ho = r oo + a ; loi = An — £ (47) 

and constrained by the following inequalities 

< [(a + /3)e] 


2[rn + hi] 
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(48) 



[e(-g + /?) - (rn - hi )] 2 

e(a + (3) — 2(rn + l\i) 


< [(a + /?)e] 


where 


ho < 


4[e(a-J-/3)-2(rn-Hn)] 


{e+*io-)'io + - 


x + P)-(rH -lll)](e + 2lQi+l 1 Q+r 1Q ) i 2 

[ e (g + g)-2(rii+i 1 i)] J 


A(\( n I /9V1 [ e (- a + /3)-(rii -ill)] 2 1 

4{l {a+PM [ e(Q+/3i _ 2 ( ril+ill )] 1 


a = 


T p- 


(3 = 


T vz 


(49) 

(50) 

(51) 


Proof : Strong stability of (46) follows if the interface treatment at x = Xi is of a dissipative nature. As we have 
already established the stability of the far-held boundary penalties we need only focus on the interface terms. 

The energy method applied to equation (46) (multiplying the discrete equations in the left and right subdomains 
by U T and V T respectively), and using the relation Q + Q T = B from the SBP rule (27), yields 


ii IlUft + ||V||» p ] 

+ 2e[||2?«U|£ | + ||2? r V|£j 


[2eUViU - aU 2 }^ + 
2loo Ut [Ut - Vi] + 
2l w (V l U) i [U i -V i \ + 
2r 00 ^ [V - U z ] + 

2r 10 (D r V) i \y i -U i ] + 


[2 eVV r V - a V 2 }} + 
2 loi Ui [( ViU), - {p r V)i\ 

2lu(ViU) i [(V l U) i -{V r V)i] 
2r 0 i Vi \{V r V)i - {ViU),} 
2r 11 (V r V) i {{V r V) i -{V l U) i ] 


(52) 


The stability of equation (52) is ensured if all the terms on the (RHS) are negative. Once again, a small portion of 
the diffusion terms ||2?jU||p and ||Z> r V||p arising at the interface can be moved to the (RHS) to help in the proof of 
stability. Defining 

Vi = Vi - ae ; _ ei_ T 


with 


V r = Vr — (3 e ; 

P = ... T 4 -1 .... 


'l ^ 

and collecting all the interface contributions into one term, (52) becomes 


A 

dt 


I|U|| 


V, 


IIVII 


2e 


P«U|| 


V, 


ID. VII 


= [2eUViU - aU 2 ]^ 







+ t 2 

eVV r V - aV 2 ] 1 

+ 

T i 

with the interface term Y,; 

= [Ii] T Mi Ii defined by 





Ui 

T 

(—a + 2/qo) 

o 

o 

8* 

+ 

o 

o 

T 

(e + l oi + ho) 

— (loi + Do) 


Ui 

T — 

J- 1 — 

Vi 


O 

O 

+ 

o 

o 

T 

(+a + 2roo) 

— (Hu + ho) 

(-e + roi +ri 0 ) 


V 

C Wh 


7t7 

+ 

o'" 1 

+ 

o 

— (Hu + ^io) 

—2(ae + hi) 

— (hi + i'ii ) 


(ViU)i 


(PrV)i 


1 

o'" 

+ 

o 

(~e + r 0 i + J*io) 

— (hi + Di) 

-2(f3e + m) 


(V r V)i 


(53) 

(54) 

(55) 


(56) 


(Again, see Appendix (C) for an estimate of the magnitude of a.) 

The first and second (RHS) terms are negative semi-definite based on single domain arguments, whereas the 
stability of the interface term T, is guaranteed if the matrix M, is negative semi-definite: 

T i = If Milt < 0. (57) 


The definiteness of a symmetric matrix is completely determined by the signs its eigenvalues, specifically if the eigen- 
values A i of the matrix M* satisfy Xi < 0 ; i = 1,4, then M t is negative semi-definite. 

The stability analysis is simplified if Mi is rotated with a similarity transformation. Define the new vector Ii = Sit 
such that: 



Ui + Vi 


'110 0' 


Ui 

i- 1 

Ui-Vi 

1 

1-10 0 


Vi 

l ~V2 

(V t U)i + (V r V)i 

"72 

0 0 11 


(ViU)i 


(ViU)i - (V r V)i 


0 0 1-1 


(V r V)i 
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The rotation matrix S replaces the stability condition given in equation (57) with the following equivalent condition: 


/■/ M"l, = XjS T SM i S 1 SX i = XfMiXi < 0; 


(59) 


where 


M, = 


0 

—a + loo - r 00 

0 

e + loi- r 0 i 


—a + loo ~ roo 
2(loo + J'oo) 
e + l io — Do 
Zoi + Zio + 1“01 + Do 


0 

e + l io — Do 
-((a + /3)e) 

-a + /3)e + In - rn 


e + ioi - r oi 
Zoi + Zio + r 0 i + Do 
(-a + 0)e + lu -r n 
— ((a + (3)e) + 2(Zh + rn) 


(60) 


Every submatrix within M t must be negative semi-definite to ensure the semi-definiteness of M$. Specifically, the 
zero at the (1,1) element of Mi, forces the first row (and column) to be entirely zero. Thus, the conservation conditions 
loo = roo + a ; and lo i = roi — e , (equation 47) are necessary conditions for 'P-norm stability of the interface. Upon 
substitution, the resultant matrix becomes 


Mi = 


0 


0 


0 


0 

0 —2a + 4Zoo e + Zio — Do e + 21q\ + ho + rio 

0 e + Zio-rio ~{(a + (3)e) {-a + /3)e + In - r n 

0 e + 21 qi + ho + rio ( — ex + 0)e + In — rn — ((q + /3)e) + 2(in + rn) 


The necessary algebra is greatly simplified if a change of variables is introduced. Defining 

t\ — (e + loi + Zio) ; Z 2 = (Zoi + Do) 
rn = e[(— a + 3/3) / 4 — Si(a + (3)\ ; In = e[(3a - /3)/4 - s 2 (a + /?)] 

and substituting these new variables into Mi yields 


(61) 


(62) 


Mi = 


0 0 0 
0 — 2a + 4Zoo t\ — f 2 

0 t\ - f 2 — [(a + /3)e] 

0 *i + * 2 [(a + /3)e][si - s 2 ] 


0 

ti + 1 2 

[(a + /3)e][si - s 2 ] 
— 2[(a + /3)e][si + s 2 ] 


Sylvester’s Theorem can be used to rotate M, into diagonal form without changing the signs of the eigenvalues. 
Repeated use of the theorem produces the following seven necessary conditions governing negative eigenvalues. 


— 2a + 4 /qo < 0 ; — [(cc + /3)e] ^ 0 ; — 2[(ct + /3)e][si + s 2 ] < 0 


[ 2a + 4Z 0 o] < [(a+foe] ’ t 2a + 4Z 0 o] < 2 r( a +mirf] 


< 1 


[(a+/3)e][si+s 2 ] ’ 2 [si+s 2 ] — 

_ 2o + 4/nn < - ftl+ t 2 i 2 

- 2[(a+/3)e][si+s 2 ] (a+(3)e\ 1- ] 

Recalling that the variables a , /? and e are positive quantities, these seven constraints can be reduced to just three 
independent inequalities. Thus, the stability of the method is governed by the following three conditions 


0 < [si + s 2 ] 


0 < 1 - 


Zoo < 


[S1-S21 2 

2[si+S2] 

[ti+^2] 2 


rj. J. I [ S 1 - s 2 ][*l+ t 2 ] i 2 

[tl-t 2 +— ( 31+32 ) I 


(63) 


2 8[(a+/3)e][ Sl + S2 ] 4 (a + /j)e[ 1 - ] 

Figure 3 shows the admissible parameter space for the variables si and s 2 , as determined from the first two inequalities. 
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Figure 3: Stability plane for parameters Si and S 2 - 


Written in terms of the original penalty parameters , equation (63) becomes 

[2(rii + Zn)] < [(a + /3)e] 

s [<«+»] 

r — i 7 I [e( — a + /3) — (rn -ln)](e+2l 01 +l 10 +r 10 ) 2 

; < a _ (e+2l 01 +l 1 o+rio) 2 _ 1 e+ho-ri 0 + [e(a+/3) _ 2(ril+lll)] -} 

00 ~ 2 4[e(a+/3) — 2(r 11 +Z 11 )] FPj±|E^II 

[e(c+/3)-2(r 11 +! 11 )] t 

the desired result for the eight parameters that ensures negative semi-definiteness of Tj in equation (55). 
Finally, integrating equation (55) in time leads to 

\\U{x u t )\\ 2 Vi + \\V(x r ,T)f Vr + / 0 T e-^- T ){ 2e[||D,U||| i + ||2> r V||| r ] + BC7 d - T<}dt 

= e+vT {\\f(xi)\\ 2 Vl + ||/(av)||p r + Jo e-^ [Data - PosDef d ]dt} 


(64) 


(65) 


Comparing equation (65) with the definition of semi-discrete strong stability [Definition (2)] leads immediately to the 
desired result. [Consult the single domain case, equations (37) - (44), for the intermediate steps in developing equation 
(65)]. □ 

Remark. The stability constraints (48), (49) relate the acceptable values of the parameters Zn and rn, relative to 
the background dissipation inherent in the discretization brought about by the parameters a, (3 and e. 

Remark. Constraint (50) describes the acceptable values of the parameter Zoo, hr terms of Zoi, Zio, Zn, j*io, and 
rn. The variables Zoi, Zio, Zn, rio, and rn are subject to no explicit stability constraints, but seriously influence the 
acceptable values of the parameter Zoo via constraint (48). Stable values of the parameter Zoo must be strictly less than 
a/2 and become progressively smaller for inappropriate choices of the other five parameters. 

Remark. Setting the parameters to the values Zn = rn = 0 removes from consideration the stability constraints 
(48), (49), and greatly simplifies constraint (50). This is the original Baumann-Oden formulation. 


3.4.3 Pointwise Stability 

Equations (44) and (65) describe the temporal evolution of the single- and multi-domain numerical solutions, thereby 
establishing strong stability of each respective formulation. Groups of numerical solution terms found on the (LHS) 
of each equation, are dominated by positive and bounded data found on the (RHS). Three distinct groups of terms 
appear on the (LHS). The groups include the solution terms ||U||p and ||V||p , the derivative terms ||X>jU||p and 
||XVV||'p , an d the boundary /interface terms BC d and Tj. The solution and derivative terms are weakly bounded in 
P-norms for finite times, while the boundary/interface terms are locally pointwise bounded. 

Note that equations (44) and (65) do not establish the pointwise stability of the solution at interior gridpoints. To 
establish interior pointwise boundedness, a discrete Sobolev inequality is required. Appendix (C) provides a proof of 
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the discrete Sobolev inequality, 

m 2 < Cl \\v\\ 2 v + c 2 \\V\5\\ 2 v 

derived specifically for SBP operators. This inequality, combined with equations (44), or (65), provides the link 
between P-boundedness of U and 2?U, and pointwise boundedness of the solution over the entire domain, and lead to 
the following theorem. 

Theorem 3 For any grid function U on (0 < Xj < l)',j = 1 ,n, and a consistent derivative operator V of rank n— 1, 
then V -boundedness of ||U||p and ||2?U||p, implies a pointwise estimate of the form 

||U||L ^ ( C mpm + Z _1 )||U|| 2 + Cmpm^WV U|| 2 (66) 

Proof : See Appendix D for the proof of Theorem 3. Extension to two domains is straightforward. 

In addition to providing a stronger measure of stability, a pointwise stability estimate can be a necessary condition 
used in proofs of global accuracy. Indeed, a pointwise estimate is the critical condition used in the work of Svard and 
Nordstrom [38], and in Section 4 of this work, to establish a sharp estimate of global solution accuracy. 


3.4.4 Primal Conservation 


We are interested in numerical solutions to the 1-D equation in divergence form 

Ut + fx = 0, |a;| < l,t > 0 ; / = au - eu x 


(67) 


where / involves convective and diffusive terms arising from for example the nonlinear Burgers’ equation. Assume that 
the discrete solution converges to a function u(x,t). We then ask the question whether this convergent limit function 
u(x,t) is a weak solution to equation (67) and satisfies 

[ <&(x)u{x)\l dx + f $(t)f(t)\t{dt = f f [$ t (x,t)u(x,t) + $ x (x,t)f(u(x,t))]dxdt (68) 

J- 1 J o J o J - 1 

for any smooth function <I>(x,t). 

The definition of conservation for the diffusive terms is somewhat ambiguous. (The original Lax-Wendroff Theorem 
was presented for the hyperbolic case, not the parabolic case.) For example, using the relation 


J $ x (x,t)u x (x,t)dx = J $ xx (x,t)u(x,t)dx - § x (x,t)u{x,t)\f_\ 


equation (68) could be further manipulated into the form 


I <&(x)u(x)\o dx+ = J J [(<&t{x,t) + a$ x (x,t) + e$ xx ) u(x,t)\dxdt 


(69) 


(70) 


A variant of the classical Lax-Wendroff Theorem is now presented in the context of the SBP formulation defined 
in equation (46). Special attention is given to the interface conditions, and the diffusive terms. 


Theorem 4 Two conditions are necessary if approximation (46) is to be conservative (in the V norm) across the 
interface. They are 

loo = ?’oo + a, l oi = ?’oi — e, (71) 


Proof : Multiplying (46) with and A>(P r . (neglecting the farfield boundary terms) leads to 


t 

+ 

&JV r Vt 










+ ai^ViVi U 

+ 

$JV r V r V] 










e[ 

+ 

A>JV r V r V r V] 

= 

loo 

[Ui - 

-Vi\ 

+ 

loi 

i’l 

mu)i - 

- (p r V)i 




+ 

ho 

(*w- 

~Vi\ 

+ 

hi 

($'] 

UlPMi- 

- C VrV)i 




+ 

roo 

Wi - 

-Ui] 

+ 

?’01 

( I’i 

[(PrV)i ■ 

- (P(t/)i 




+ 

no 

(^O JK: - 

-Ui] 

+ 

ni 

($'] 

h[V>rV)i - 

- (ViUfi 
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( 72 ) 



where <f>, and (<J >/ ) . are the values of the test function and its derivative at the interface location x,;. Both are continuous 
across the interface. 


Rearranging the time terms yields 

QfPlUt + QjVrVt = + *TP r V) - P lf t Vl U + ($ r )fp r V] 

Making use of the SBP rules (26 - 27) on the advection terms yields 

QjVlDiU + &?V r V r V = -i&tfViU - (&r) T VrV 

+ [mti + mt 1 


(73) 


(74) 


Again, note the ambiguity in the meaning of conservation for the diffusive terms. Using the SBP rules (26 - 27) once 
on the diffusion terms yields, 


GfVlViVt U + ^V r V r V r V 


-WfpfoU) - {® r ) T v r v r v) 

[*«(ZW)]|-1 + [*r(PrU )] |£ 


while twice yields 

^VlVW + ^V r V r V r V 


(V/fPlU) 

mviU) -WM-1 


+ WfVrV) 

+ [3> r {VrU)-{& r )V\\£ 


(75) 


(76) 


Inserting (73), (74) and (75) into (72), yields the condition 


f t ^T ViXJ + ^T VrY) 

+ [^ l (a~eV l )U] |_r + [$ r (a-eV r )V] | +1 


mft-PlU + ($r)JV r V] 
{& l ) T Vi(al l -eV l )\J + {& r ) T V r {al r 
&i{loo — r oo — a){Ui — Vi) 

$<(ioi -roi + e)(DU l -VV l ) 
mho -r 10 ){Ui-Vi) 
mhi -rii)(VUi-VVi) 


eV r )V 


(77) 


Integrating equation (77) with respect to time, applying the conservation conditions given by equation (71) and the 
additional conditions ho = ?’io and hi = bi produces the equation 

(*TP l V + *?P r V)\'S 

+ / 0 T {[4> i (a-cA)C/]|_i + m-eV r )V}\+ l }dt = ^[(^^U + ($ r )?> r V]di 

+ Jo fii&ifViiah - eP/)U + ( <V r fv r (al r - eV r )\}dt 

(78) 

which matches the weak solution given by equation (68) term by term, in the limit of infinite spatial resolution. □ 
Inserting (73), (74) and (76) into (72), yields the condition 


^fVlV+^VrV) 
+ [^ l (a-eV l )U] |-r + [$ r (a-eV r )V] | +1 


mfiPlU + mVrV] 

{<Vi) T Vi{ah) U + ( <V r ) T V r {aI r )V 

ei&hfViU + e($" r ) T U r V - e[$'C7 |_i + $' r U | +1 ] (79) 

$i(ioo -roo - a)(Ui-V) 

$i(/oi -r 0 i + e)(DUi - VVj) 
mSio - r 10 - e)(Ui-V) 

(&)i(hi - ru)(T>Ui-VVi) 


Integrating equation (79) with respect to time, applying the conservation conditions given by equation (71) and the 
additional conditions Zio = rio + e and hi = hi, we see that expression matches the weak form of the advection 
diffusion equation in the limit of infinite spatial resolution. 

Remark. Note that the auxiliary constraints in equations (78) and (79) depend on the definition chosen for 
conservation. Still another definition (perhaps a linear combination of the two) would have produced a third set of 
constraints. 


16 



3.4.5 A ”flux-form” discretization technique 


A finite difference approximation of the coupled problem (168) motivated by the Local Discontinuous Galerkin FEM 
proposed by Cockburn and Shu and referred to as “LDG” [14] is 


ViU t + aViViV — eViT>i<j) 

= L m [Ui-Vi] 

ei_ 

+ 

LoMi 

- Wi. 

eVi (</» - ViU) 

= L w [Ui-Vi] 

ei_ 

+ 

LoiMi 

- Wi. 

U(s,0) 

= 0. 





V r V t 4“ aP r UfV — eP f T2 r 'ip 

= Roo[Vi-Ui ] 

e i+ 

+ 

RodWi 

- (4>)i. 

eV r ( - V r \) 

= R 10 [Vi-Ui] 

e i+ 

+ 

RnWi 


V(*,0) 

= 0. 






Theorem 5 The approximation (80) of the problem (168) is strongly stable if the eight parameters are related by the 
two equalities 


Loo — -Roo + « ; L oi — Rqi — e 


and constrained by the following inequalities 


2[_Rii+Lii] < [(a + P)e] 


(81) 


(82) 


[e (— a + (3) — (i?n — in)] 
e(a + (3) — 2(l?ii + in) 


< [(a + /3)e] 


(83) 


^oo ^ 


< # - 


(e+2Loi+I/io+-Rio) 

4le(a+/3)-2(flu+L 11 )J 


where (as with the primal form) 


r_,r r, , [e(-c+/3)-(K 1 i-in)](€ + 2I, 01 +L 10 + B 10 ) ,2 

{e + L 10 -Kl0+ [e(c+/3)-2(R 11+ Z. 11 )] 1 

A t I 0\ . [e(-c+/3)-(-Rn -in )] 2 

4{(a+/^)e [« Cci+ ^_ 2 (k 11+£i1i) ] 1 


a = 


/3 = 


(84) 


Furthermore, the stability constraints for the LDG scheme are identical to those of the primal scheme although stability 
is proven in different variables. 

Proof : Strong stability of (80) follows if the interface treatment at x = Xi is of a dissipative nature. As we have 
already established the stability of the farfield boundary penalties we need only focus on the interface terms. 

The energy method applied to equation (80) is derived by multiplying the first and second discrete equations in 
the left subdomain by U T and <p T , respectively, and the first and second discrete equations in the right subdomain by 
V T and i p T , respectively, and using the relation Q + Q T = B from the SBP rule (27), yields 


a iiuft + iivii’j 

+ HM 2 v, + Mv r ] 


[2 eU(f> - aU 2 ] ^ + 
2 L 00 Ui [Ui - Vj] + 

2L w (<t>) i [U i -V i } + 
2R 0 o Vi [Vi - Ui] + 
2RioMi[Vi-Ui] + 


[2 eVxf - a V 2 ],) + 
2L 01 Ui {(f)i - Mi] 
2L 11 (4>)i[(<t>) i - 
2 Roi V [Wi - (0)i] 
2 Ri 1 (tp)i[(ip) i - (</>)*] 


(85) 


The stability of equation (52) is ensured if all the terms on the (RHS) are negative. 
Collecting all the interface contributions into one term, (85) becomes 


d 

dt 


l|u|£, 


IV 


2e 


Mp 


Vi 



with the interface term T, = 7) T M, 7) defined by 
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[2 eU(j> - aC7 2 ]_ 1 
[2 eVfi - a V 2 ] 1 + Ti 


(86) 



(87) 


T, = 


Ui 

T 

(—a + 2 Loo) 

— (Loo + Roo ) 

(e + L01 + L10) 

— (Lqi + Rio) 


' Ui n 

Vi 


-(Too + Roo) 

(+a + 2 R 00 ) 

— (Roi + bio) 

(— e + Roi + Rio) 


Vi 



(e + L01 + L10) 

—(Roi + Fro ) 

— 2 (ae + Ln) 

— (T11 + Ru) 


(<t>\ 

(VO* . 


— (Lqi + Rio) 

(— e + Roi + Rio ) 

— (T11 + -R11) 

— 2 (( 3 e + Ru) 


_(V 0 ; . 


Comparing equation (87) with (56), we see that the stability matrix M, for the LDG formulation is identical to 
that obtained with the primal formulation. Applying similarity rotations produces an LDG conservation condition 

Too = Roo + a ; Lqi = Roi — e, (88) 

Manipulating the matrix M. j, (exactly as in the primal formulation), yields the stability inequalities 

0 < 

0 < 

Loo < 


Figure 3 again shows the admissible parameter space for the LDG variables Si and S 2 . Thus, the LDG and primal 
stability conditions are identical. Transforming equation (89) back to the original variables L 00 — R 11 yields the desired 
result for the eight parameters that ensures negative semi-definiteness of in equation (86). 

Finally, integrating equation (86) in time leads to 

\\U(xi,T)\\ 2 Vi + \\V(x r ,T)\\% r + J^e-^-U{2e[\m% + IMk] + BC ^ ~ (QQ) 

= e+r < T {\\f{xi)\\ 2 Vi + \\f(x r )\\ 2 Vr + Jq e-^ [Data - PosDef d ]dt} 

Comparing equation (90) with the definition of semi-discrete strong stability [Definition (2)] leads immediately to the 
desired result. □ 

Remark. Although the stability boundaries for the eight interface parameters are identical for the primal and LDG 
schemes, as are the norms used in both cases, the two methods are stable in different diffusive variables. 

Remark. As with the primal formulation, the stability constraints (82,83) strongly constrain the acceptable values of 
the parameters Ln and Ru, as compared to the background dissipation inherent in the discretization. The constraint 
(84) describes the acceptable values of the parameter Loo, in terms of L 01 , L 10 , Lu, Rio, and Ru. The variables L 01 , 
L 10 , L 11 , i?io, and Ru are subject to no explicit stability constraints, but seriously influence the acceptable values of 
the parameter Loo via- constraint (84). 

Remark. As with the primal formulation, a proof of pointwise stability is achieved using a Sobolev inequality. We 
state without proof, the following theorem. 


[Si + S 2 \ 

l _ [Si-S2l 2 
1 2[S;l+S 2 ] 


[Ti+T 2 ] 2 

8[(a+^)e][S!+S 2 ] 


r-T- ^ , [Sl-S 2 ][Ti+T 2 h2 
[h-l2+ 2 ( Si +S 2 ) — I 

4(a+/3)e[l— ] 


(89) 


Theorem 6 For any grid function U on (0 < Xj < l)',j = 1 ,n, and a consistent derivative operator V of rank n— 1, 
then V -boundedness o/||U||p and ||4>||p, implies a pointwise estimate of the form 

||U||L < (Cmpm + Z-^HUH 2 + C mpm - 1 ||$|| 2 (91) 

Proof : See Appendix D for the proof of theorem 6. Extension to two domains is straightforward. 


3.4.6 LDG Conservation 

Again, we are interested in accurate numerical solutions to 1-D divergence equations of a form consistent with that of 
equation (67). We then ask the question whether the convergent limit function u(x,t) is a weak solution to equation 
(67) and satisfies equation (68). The convective terms at the interface in the LDG formulation are penalized the same 
as in the primal formulation. It is reasonable to expect the necessary conservation conditions to be the same in the 
two formulations. The diffusive interface terms are very different between the two formulations. We now prove the 
conservation theorem in the LDG formulations. 
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Theorem 7 Two conditions are necessary if approximation (80) is to be conservative (in the V norm) across the 
interface. They are 

Loo = Roo + a, Lqi = i? 0 i — e, (92) 


Proof : Multiplying (80) with <&JVi and V r , 

[using the SBP rules (26 - 27)] leads to 


(neglecting the physical boundary terms) and rearranging terms 


+ *TVrV r l)> = -(VfVK)) - (Kf'Prlp) 

+ + i^mt 1 


(93) 


is the value of the test function at the interface location ay, and it is continuous across the interface. 
Manipulating (73), (74) and (93) yields the condition 


|($^U + *?P r V) 


imf-PlV + ($r)JVrV] 

(<f>'i) T Vi{(aU - e</>) - eU] + ($' r ) T R r [(aV - eif) - eV] 
$ z [(-af/ + €(/>) - eU] |_i + $ r [(-aV + ei)>) - eV] |+i 
^i(^oo — Roo — a + e + Lio — Rio)(Ri — Vf) 

*&i{Loi — Roi + e + Ln — Rii){4>i — ifi) 


(94) 


Integrating equation (94) with respect to time, applying the conservation conditions given by equation (92) and 
the additional conditions L io = -Rio — e and In = flu produces the equation 

(* 1 ^ u+^v)u 

+ / 0 T {[4> i («^eR i )C/]|-i + ['5>r(a-eV r )V}\ +1 }dt = /J[($,)fPj U + (<f> r )fV r V}dt 

+ JoU&if'PiiaZi - U + (<S>'r) T Vr(aI r - eV r )\}dt 

(95) 

which matches the weak solution given by equation (68) term by term, in the limit of infinite spatial resolution. □ 

Remark. Only one set of constraint equations results from the conservation condition in the flux formulation. 
Further, the constraints are different than those obtained by either definition of conservation in the primal form. 


3.5 Relating the Primal and Flux Formulations 

One can imagine the possible existence of other formulations in addition to the primal and flux collocation formulations. 
On the other hand, one might also imagine that both are equivalent (as shown in FEM by Arnold et al [3]) with an 
appropriate change of interface variables. We begin by showing that the LDG formulation can be implemented in 
terms of the primal scheme if specific values of the interface parameters are used. 

Theorem 8 The LDG interface parameters Loo - Rn are related to the primal interface parameters loo ~ hi by the 
following relations: 


^00 

= 

Loo 

+ LoiCi 

+ 

LlO 

a 

+ 

Luci 

a 

hi 

= 

0 

+ L 01 C 2 

+ 

0 

+ 

L 11 C 2 

cx. 

ho 

= 

0 

+ 0 

- 

Lio 

- 

LuCi 

hi 

= 

0 

+ 0 

+ 

0 

— 

L 11 C 2 

Do 

= 

O 

o 

+ R 01 C 1 

- 

-Rio 

0 

- 

fluci 

(3 

r 01 

= 

0 

+ R 01 C 2 

+ 

0 

- 

-Rl 1 C2 

0 

no 

= 

0 

+ 0 

- 

Rio 

- 

R 11 C 1 

rn 

= 

0 

+ 0 

+ 

0 

— 

R 11 C 2 


Cl 


C2 = 


1 




i-h 



(96) 


(97) 
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Proof : The LDG two-domain formulation is given in equation (80). Solving equation (80) for <j> and -0 yields the 


expressions 


= ViU + \V^{L w [Ui - Vi] 

1 -TJ — 1 


e i + L n [(4>)i - WOJ e i- } 

= VrV + IPrHRloWi-Ui] e i+ + Rn [{4>)i — e i+} 

while, solving equation (98) for the interface values yields the expressions 


= (Vi\J)\i + ±{L w [Ui-Vi\ + LuK^i-Wi]} 

1>i = e i+ T V- = (PrV)\i + ^{Rio[Vi-Ui] + Rum i -(^ i }} 

with the difference being 

4>i -l>i = (X>,U)U - (PrV)\i + + + |(^ + ^)[0i-^i] 

Rewriting equation (100) in terms of the difference <f>i — 'i/’i yields 


k - if>i = Cl [Ui-Vi\ + c 2 [(: D,U)|< - (2? r V)|i] 

1 


ci = 


1 Lip | Rip ^ 


wi^F+^F) 


C2 = 




(98) 


(99) 


( 100 ) 


( 101 ) 


Substituting equations 
expression 

ViUt + aViViXJ 


(98), (99) and (101) back into the original LDG formulation [equation (80)] yields the 


P j ■ V / + a'P r T) r '\F 


= eTiViViU 

+ {L 0 o[t/. - V] + L 0 i( Cl [E/i - VS] + C2[(2?,U)|j - (2> r V)|i])} 
+ {L w [Ui - VJ] + Ln(a[Ui - Vi] + c 2 [(D;U)| i - (2? r V)|i])} 

= eV r T> r T > r'V 

+ {R 0 o[Vi - Ui\ + RoMVi - Ui\ + C2[(2? r V)|j - (£> r U)|i])} 
+ {R w [Vi - £/i] + Rn(ci[Vi - U t ] + c 2 [(D r V)|i - (2? r U)|i])| 


Rl Gi 

ViViVfW 


(102) 


l r 

V r 'D r 'Pp l e 


J1 + 

4 


Simplifying equation (102) using the relations 


ViViVfW 

Vr'Dr'Pp'e 


!+ 


= /VG + 

= -Vj ej - 


— ei 

(3 t ' I + 


yields 


ViU t + aViViXJ - eViViV , U = 

+ 

+ 

+ 

P r V t + — eR r D r 7b r V = 

+ 

+ 

+ 


-boo 

Rio 


-boiCl 

Rll^l 

a 

+ -boiC 2 

-^llC2 


0 

0 + 

<JL 

0 + 0 
— Lio — -LnCi 
0 + 0 
0 — LuC2 

Roo + Roici 

-Rio RllCj 

P P 

0 + -R 01 C 2 

n Rj 1 a 

u 0 
0 + 0 

—Rio — R 11 C 1 
0 + 0 
0 — Ri 1 c 2 


[{/* - ^] Xi e ; _ 

[(x>iU)U - (Z> r v)|i] h ei_ 
[Di-V] X>; T e;_ 
[(X> Z U)U - (2? r V)|j] Zfei_ 

Wi - Ui ] lr e i+ 


(103) 


] [(2? r V)|i - (2? r U)|i] l r 

] [V - Ui] vje 

} [(2? r V)|i - (2? r U)|i] Vje, 


- 1 + 


*+ 


Comparing equation (103) 

Remark. Note that the 
the new parameter 


with the original primal scheme given in equation (46) gives the desired result. □ 
primal coefficients are not constants, but rather change with grid resolution. For example, 


^00 — -boo + ioici + — - 
a 

20 


-bllCl 



given in equation (96) is inversely proportional to the grid spacing Ax since a oc Ax, and becomes larger as the grid 
is refined. 

It is reasonable to ask whether the stability constraints (48) and (49) for the primal method are preserved for 
methods satisfying the LDG stability constraints (82) and (83). 

Theorem 9 Stability (in the V-norm) in the LDG variables is a sufficient condition for stability (in the V-norm) in 
the primal variables. 


Proof : We begin by assuming that the LDG variables Loo — > Rn are chosen such that the conservation conditions 
given by equation (88) are satisfied 

Loo = Roo + a, L oi = Roi — e, 
that values of the parameters Si and S 2 satisfy the constraint equations (89) 


0 < [S 1+ S 2 ] ; < 1 

and 

* ~ . r.S i-S 2 lfT 1 +T 2 l 1 2 

2(Si+S 2 ) 1 

_ [S 1 -.S 2] 2 
2[Si+S 2 ] 

Substitution of the transformation relations given by equation (96) into the primal conservations conditions loo = 
roo + « ; and /qi = Lji ~ e , yields the conservation equations based on the LDG variables; 


L 


00 


< 


[Ti + T 2 \ 


- ± 2 + 


8[( a + /3)e][5'i + S 2 ] 4(a + /?)e[l 


^00 — r 00 — a , 
loi ~ roi + e 


0 

0 


(Loo ~ Roo — a) + 


4:(aRio+0Lio)(Loi — Roi+e) 

€[(a-/3)2+4(a+/3)(a5i+/3S 2 )] 
4(a^)(L 0 i — flpi+e) 

[( a _ / 3)2 +4 ( Q , +/3 )( QSl+/ 3s 2 )] 


0 

0 


(104) 


while substitution of the transformation relations given by equation (96) into the primal and the stability constraints 
given by equation (63) or equation (64) yields the stability equations based on the LDG variables; 


and 


0 < [si + S2] 


2[a 2 (l+45 1 )-/3 2 (l+4S 2 )] 2 +64a 2 /3 2 [5 1 +S 2 ][l-fi^|i y 

[2((a-/3) 2 +4(a+/3)(aS 1 +/3S 2 ))] 2 


0 < [1 


(■S 1 -S 2) 2 1 
2(si+S2) ^ 


0 < 


32a 2 /3 2 [S 1+ S 2 ][l-f^||^] 

[a 2 (l+4Si)-/3 2 (l+45 2 )] 2 


loo < f 

-boo < f 


l^i +^2! 2 

8[(a+/3)c][si+s 2 ] 

[T!+T 2 ] 2 

8[(a+/3)e][Si+S 2 ] 


4(a+« e [l-^j^T 


rr . , [Sl-S 2 ][T 1 + T 2 ],2 
~ 1 2 +Sj] i 


(105) 


(106) 


Thus, if the original scheme written in the LDG variables satisfies the required conservation and stability constraints, 
then it also satisfies the corresponding conservation and stability constraints written in the primal variables. □ 


4 Generalized Penalties 

Equation (46) presents a new procedure, (Baumann-Oden FEM), for imposing penalties at the interface between two 
subdomains or elements. The penalties are constructed from precise combinations of “jumps” in the solution and first 
derivative across the interface. Note, however, that in general, penalties could be constructed from any smooth quantity 
that spans the interface. Thus, reasonable candidates for penalties include the solution and its first p derivatives. 
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A semidiscretization of equation (30) that includes penalties on the solution and the first two derivatives written 
for the two domain problem (in BO form) is 


ViUt + aViDiU = 


U(®, 0) 
Vr'Vt o'Pr 7-V V 


V (x, 0) = 


(PiViVi U 

{loo[Ui-Vi\ + l 0 i[(VlU) 
{l w [Ui - Vi\ + l n [{VjU) 
{lw[Ui-Vi\+hi [iV]U) 
{a_ 1 (L? 1 (U))e_ 1 } 

0. 

eVrVrVrV 

{roo[Vi - Ui } 

{riolVi-Ui] 

{?’20 \Yi - Ui] 

W+i (L+i (V))e +1 } 

0. 


- {VlV)i] + *02 [(VfU)i - (PrV)i]} [h) T 

- <PlV)i] + iMCDfU), - {T>rV)i\} m 1 

- ( 2 Wi] + h 2 [(Vfu^ [ 2 ??]- 




(Piu)i\ 

(PlU)i\ 


■ro 2 [(PWi-pfUWpr 


roi[(.'D 1 r V) i 

rn[(^V)i 

T 2 i[{VlV) i - (VjU) i )+r 22 [(V* r V)i - (VfU)i}} [Vl 


■ri 2 [(ViV)i-(VfU)i]}[V r 


-1+ 

e i+ 




(107) 


4.1 Stability 

Theorem 10 The approximation (107) of the problem (30) is strongly stable if the eighteen parameters are related by 
the six equalities 

loo = J'oo + a ; loi = r 0 i — e 

^02 = I02 ; ^20 = ho ; ^21 = hi ; t 22 = l 22 


and constrained by the following inequalities 

2[rn + Zn] < {(a + j3)e] 


(108) 


(109) 


[e(-a + P) - (rn -/n)]~ 
e(a + (3) — 2(?’n + /n) 


< l(o: + P)e] 


( 110 ) 


with 


and 


once again 


< - 


122 


loo < 


(+2r2i+ri2+^i2) 
4[e(o!+/3) — 2(rii+Jn)J 


[ti+t 2 ] 2 

8l(a+/3)e][si+s 2 ] 


111 _ N, [«(-a + /3)-(»'ll -<ll)l(+ 2 »'21+''12+ ! 12) 1 2 

l(h 2 -ri 2 )+ [e ( Q+g )_ 2 ( r . 11+ill ) 1 } 


mw~ i o\ _i Lei— a-t-pj-im — n i JJ 1 
4{[(a+/3)e] [ e(c<+/ 3 ) _ 2 ( ril +ill) ] 1 

r. . i I s ! -« 2][*1 +* 2 ] l 2 

[«1~*2+ 2 (s 1+ s 2 ) 

4(a+/3)e[l-^f^T 


r „ T , . ^- t2 + ui ^&r 2 ^ -»2 + t %? 1 ) &r 2) 1 1 2 

{ 1+ ™™ + ^ } 

~ ^ 1,1 
-4{i22 + 


2L S 1 + s 2 J 


t\ = (e + ^oi+^io) ; O = (^oi+Oo) 
r\i = e[(-a + 3/3)/4-si(a + / 3)] ; bi = e[(3a - /?)/4 - s 2 {a + /3)], 

24 = {r 20 + l 02 ) ; Vi = (^ 12 +?’ 2 i) ; J /2 = {ri 2 + r 2 i) 


a = 


/3 = 


( 111 ) 


( 112 ) 


(113) 


(114) 


Proof : Strong stability of (107) follows if the interface treatment at x = Xi is of a dissipative nature. As we have 
already established the stability of the farfield boundary penalties we need only focus on the interface terms. 
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The energy method applied to equation (107) by multiplying the discrete equations in the left and right subdomains 
by U T and V T respectively, yield 

U T Vi£lJ + V T V r £V = -a[U T Q l U + V T Q r V]+e[ IF QfoU + V T Q r V r V] 

+ Ui {Zoo [Ui - V] + l 0 i[(ViU)i - (V r V)i) + l 0 2[(V?U)i - (V 2 V) t ]{ 

+ (V}U)i {l 10 [Ui - Vi\ + Jn[(ZW)< - (®rV)i] + h 2 {(VfU)i - (V 2 r V)i}} 

+ (V?U)i {ho [Ui - v;] + hi mu)i - (V r V)i] + h 2 [(V 2 U) i - (V 2 V)i}} (115) 

+ V {r 00 [Vi - Ui] + r 0 i[(V r V)i - (2W)i] + r 02 l{V 2 r V)i - (V 2 U),}} 

+ (VlV)i { n 0 [Vi - Ui] + r u [{V r V)i - (2W)i] + r 12 [{V 2 V)i - (VfU).]} 

+ C V 2 r V)i {r 20 [V - Ui] + r 21 [{D r V)i - (2W)i] + r 22 [(V 2 V)i - (V 2 U).]} 

Using the relation Q + Q T = B from the SBP rule (27), borrowing a small portion of the ||X?fU|| p and ||XVV||p 
terms arising at the interface, and collecting all the interface contributions into one term, equation (115) becomes 

Tt [||U|£ ( + ||V|| = J + 2e [||2W|& + ||2? r V||| p ] = [2 eUV t U - aU 2 ]_ ± (n0) 

+ [2eVD r V - a V 2 ] 1 + T 4 

with the interface term T , = [Ji] T Mi Ji defined by 

Ji = [ Ui, Vi, {T>iU) i , (T> r V) v {V 2 U) v (V 2 V) i f (117) 


(— a + 21 qo) ~(Zoo+?’oo) (e + Zoi + ^io) — (Zoi + **io) (Z02 + Z20) — (^02+^20) 

-(loo+roo) (+a + 2 r 00 ) — (nn+Zio) (-e + r 0 i + ri 0 ) -(£20+^02) (r-o2+r 2 o) 

(e + Z01 + Uo) —{i'01 + ho) 2 (Zn — ae) — (Z11+JT1) (Z12 + Z21) — {I12+X21) 

-(loi+no) (-e + roi+rio) -(Zii+rn) 2(m - /3e) ~(hi+n 2 ) (r 21 +r 12 ) 

(Z02+Z20) —(ho + 1'02) (Z12+Z21) — (^21+^12) 2l 2 2 —{ 122 +^ 22 ) 

~ (Zo 2 +?’ 2 o) (ro 2 + r 2 o) ~ (Zl 2 +?’ 2 l) (j’ 21 +^ 12 ) — (^ 22 +^ 22 ) 2 7*22 

(H8) 

Stability of the interface is guaranteed if the matrix M t is negative semi-definite which yields an interface term of 
the form Yj = [J)] T Mj < 0. The definiteness of a symmetric matrix is governed by the sign of its eigenvalues. 
Thus, Sylvester’s Theorem is used to systematically rotate Mj into diagonal form while maintaining the signs of its 
eigenvalues. 

Applying similarity rotations, yields the generalized conservation conditions 

ho — r oo + « ; hi = i'oi — e 

ro2 = I02 ; ^20 = ho ; P21 = hi ; t 2 2 = I22 

and changing variables 

ti = (e + Zoi+Zio) ; t 2 = (Zoi+^ 10 ) 

ni = e[(-a + 3 / 3 )/ 4 - si (<* + /?)] ; lu = e[( 3 a - /?)/4 - s 2 (a + /?)], ( 119 ) 

xi = (r 20 + h 2 ) ; 2 /i = (Z 12 +?’ 2 i) ; 2/2 = ir 12 + r 2 i) 

yields the reduced matrix 

0 0 0 0 0 

2a + 4Z 0 o ti — t 2 ti 1 2 0 2 xi 

ti-t 2 -((a + /3)e) (a + P)e(sl - s2) 0 j/i — 2/2 

h+t 2 (a + /3)e(sl — s2) — 2 (a + /3)e(sl + s2) 0 yi + y 2 

0 0 0 0 0 

2 X\ 2/i - 2/2 2/i + 2/2 0 4 Z 22 
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Repeated application of Sylvester’s Theorem produces the stability inequalities 


0 < 
0 < 


[si + s 2 ] 

1 _ [si — s 2 ] 2 
2[si+S2j 


1 22 ^ 


< 


[yi+j/ 2] 2 

8[(a+/3)e][si+s 2 ] 


lyi-y2+ lsi Zf + X) y2] f 

4(a+/?Ml-Jfcfg£] 


^00 < 


[^ 1 +^ 2 ] 

8[(a+/3)e][si+s 2 ] 




( 121 ) 


{2 Xl + 


(O ±1 ^y 1 + y, X , + <S1 + ( ; 1 7 ^:^ I 

2(a+ ^ (si+s2) 4 ( . +Wg[1 -fe^gr 


, ,2 \v 1 t; | ~1 - e 2][«i+«2] i 2 

tti 1 [»l+» 2] 2 , M1 “ 2 2(31 + 52) 

4p 22 + 8[(q , + /9H[ „ 1+S2] + [3,-3 2 ]2 

4(a+p;e[l— 2 [ 31+ , 2 | 1 


2 

} 


Changing back to the original variables [via equation (119) ] yields the desired result. □ 

Remark : A general interface procedure would include derivative terms up to p th - order, and would result in the 
interface parameters Zoo — ► l PP , and ?’oo — > r pp . A stability analysis could be used to identify conservation conditions, 
as well as the parameter space that guarantees L 2 and pointwise stability. 

Remark. A general interface procedure could be constructed based on the “flux-form” of the equation. 


5 Formal Accuracy 

We derive error estimates for both semi-discrete discretization methods. We begin with the primal form representation 
[equation (46)], and conclude with the flux form representation [equation (80)]. Extensive analysis of the accuracy of 
FEM is available in the literature (see for example the work of Arnold et al. [3]). Thus, the focus of attention herein 
is high-order finite difference methods. 

Define the vectors u = [u(x 0 , t) , • • • , u(x n , t)] T , u x = [u x (x 0 , t) , • • • , u x (x n , t)] T , and u xx = [u xx (x 0 ,*),---, u xx (. x n , f )] T , 
to be the projected values of the exact solution and derivatives, defined in the left domain at gridpoints x. Similarly, 
define the vectors v, v x and v xx in the right domain. With these definitions, the continuous problem (30) can be 
represented in vector nomenclature as 


U t 

+ 

CL\1 X t^-xx 

= 

o"" 4 

O 

IT 

1 

Vi] 

vr l h e ; 

+ 

Z 01 [u x 

- V x ]i 

ei 




+ 

ho[ui — 

Vi] 


+ 

hi[u x 

- V x \i 

V^Vfe^ 



u(x, 0) 

= 

0. 







V t 

+ 

av x - ew xx 

= 

roo[vi - 

Ui ] 

V^lr e i+ 

+ 

roi [v x 

- U x ]i 

V-'lr e i+ 




+ 

rio[vi - 

Ui] 

V- 1 Vje i+ 

+ 

ru[v x 

- U x ]i 

"P” 1 Vje i+ 



v (z, 0) 

= 

0. 








with suitable initial and boundary data, and the source term F(x, f) = 0. Note that the interface terms [iq — Vi] and 
[u x — v x ]i are identically zero for the exact solution, but have been added to the right hand side of equation (122) to 
facilitate subsequent manipulations. 

Substituting the discrete accuracy relations u^, = T>iu + T e i; and u xx = V{Divl + T e2i [see definitions (25) and 
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(28)], into equation (122) yields 


u t + aVi\i - eViViU = 

O 

O 

■jr 

i 

V^li e ; 

+ 


- (V r v)i] 

v^h e; 

+ 

o 

■jr 

i 

c2 

Pf'Pfei- 

+ 

/n[(D r v) i 

- (2?,u)i] 


- 

oTgi; 


+ 

M(Tei;)i 

- (Telr)i] 

Pf 1 h ei 

+ 

eT 6 2 1 


+ 

ill[(Telj)i 

- (Telr)ij 

Pf 1 V? e; 

u(a:, 0) 

0 . 






v + aD r w — eD r T> r v = 

r 00 [vi - u-i ] 

V-'Ir e i+ 
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?’01 [('D r v)i 

- (2?iU)i] 

Pp 1 li ei 
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rio [vi - m } 

V~ x Vj e i+ 
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rii[(T>iu)i 

- (Prvjjj 

V^Vfei 

- 

aT e l r 


+ 

r 0 l [(T e lr)i 

- (Tel,)i] 

P r 1 l r ei 

+ 

eT e 2 r 


+ 

ril[(T e l r )i 

- (Tel,)i] 

V~ x Vj ei 

v(x,0) = 

0 . 







( 123 ) 


Next, define the semi-discrete error vectors Hi = U u and H r = V — v in the left and right domains, 
respectively. Premultiplying equation (46) by V~ l and subtracting equation (123) produces an evolution equation for 
the error vectors. Specifically, the collocation error of the coupled problem (168) written in primal form is 


+ aViSi — eViViSi 
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V^Vfe^ 
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aT e i i 
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i0l[(T e il)i — (Tel r )i] 
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£ T e 2 1 
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K 1 V? e ; _ 
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^00 [Sri _ S(j] 

V~ l l r e i+ 

+ 

r 01 [(V r E r )i - (D r Zi)i\ 

Pp 1 ^ e i+ 

■ 

'O' 10 [^ ri ^ li\ 

Pp 1 Vje i+ 

+ 

r 11 [{'D r E r ) i - {D r Ei)i} 

Pp 1 V^e i+ 

+ 

aT e l r 


- 

^0l[(T e lr)i — (Telj)i] 

Pp 1 l r e ; 

- 

eT e2r 


- 

^ll[(T e lr)i — (Teli)ij 

V~ x V T r e ; _ 

= 

0. 






(124) 


Four distinct types of discretization errors appear in equation (124). The first two truncation error vectors: T e i 
and T e2 , arise from errors associated with the approximation of the first and second derivative terms, respectively. Like 
conventional HOFDM operators, they are seldom uniformly accurate throughout the domain. Points near boundaries 
are commonly discretized less accurately than those in the interior. Furthermore, the boundary stencils used in first 
and second difference operators are frequently of different orders. 

The last two error vectors result from the interface penalty terms and have the forms: [(T e i)i — (T e i)i] V -1 1 ei 
and [(T e i)j — (T e i)i] V _1 V T ei. Note that the accuracy of these vectors is influenced both by the truncation error 
of the interface derivative operators, (T e i r )j, and by the size of the penalty scaling terms. To assess the size of the 
scaling terms, first note that the differentiation operator and its transpose, [D and V T ] scale as 0(Ax)~ 1 . Further, 
since Q is unitless, and T> = V~ 1 Q , the inverse norm operator V -1 scales as 0( Ax)~ . Thus, the penalty error vectors 
scale as 0( Ax) r ~ and 0{Ax) r ~ , respectively, where r is the local accuracy of the interface derivative operator. 

In summary, error vectors in both domains have leading order truncation terms that scale as 
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= [0( Ax rl ) ,- 

■ ■ , 0(Ax rl ) 

, 0(Ax 2p ), ■ ■ 

• , 0(Ax 2p ), 0(Ax rl ) ,• 

• • , 0{Ax rl ) 


T e2 

= [0( Aaff 2 ) 

■ • , 0(Ax r2 ) 

, 0(Ax 2p ), ■ ■ 

• , 0(Ax 2p ), 0(Ax r2 ) ,• 

■■,0( Ax r2 ) 

(T e] 
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= [0(Ax rl ~ 1 ), ■ 

■■,0( Ax rl ~ 
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•, 0 , 0(Ax rl " 1 ), ■ ■ 
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■ ■ , 0(Ax rl ~ 

2 ), 0 

•, 0 ^(Aaff 1 - 2 ),-. 

■ ■ , 0(Ax rl ~ 2 



(125) 


where the r 1 and r2 exponents denote the boundary stencil order of accuracy in the first- and second-derivative 
operators, respectively. (Here, we assume the interior stencils to be the same order of accuracy for both the first- and 
second-derivative operators.) 

Local lower order error terms do not necessarily decrease the global formal accuracy. The impact of lower-order 
terms on global discretization accuracy is a long-standing area of research that dates back to the pioneering work of 
Gustafsson [18]. There it was established that hyperbolic operators could accommodate local terms one order lower 
than design order, while still maintaining formal accuracy. Recent work of Svard and Nordstrom [38] has extended this 
result to include the impact of low-order stencils for parabolic problems. Their work on global accuracy is encapsulated 
in the following theorem. 
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Theorem 11 If (46) is a pointwise stable discretization of the continuous problem (30) for Ax < Axo, then with 
interior operators V and W satisfying discretization orders of 2p — 2 < rl,r2 < 2 p at the boundary, and interface 
penalty terms satisfying discretization orders of2p—2 < r 1 — 2 < 2p , then the global order of accuracy of approximation 
(46) is 2 p. 

Proof : See Svard and Nordstrom [38], Theorems 2.6 and 2.8, pp. 335-337 for the original theorems and proofs. 

Now consider a realistic combination of operators and study the effects of each error term on global accuracy. 
Assume that the derivative operator T> is defined as in equation (125) with rl = 2p — 1 . Furthermore, assume that 
the second derivative operator is formed as the action of two first derivative operators W. Finally, assume that the 
derivative used in the penalty (obtained from the derivative operator T>) is O) Ax 2 ^ 1 ), and is used in all penalty terms 
where necessary. In this scenario, the derivative operators T> and VD maintain the design accuracy 2 p of the method, 
as well as do the penalty terms producing errors of the form (TeijiP^ 1 ■ The remaining penalty terms producing 
errors of the form V T ei decrease the formal accuracy to 2 p — 1 , at least in principle. 

Remark. Adding penalties of the form (T el ) i V~ 1 V T ei is ill advised on a theoretical basis, although experiments 
presented later show that these estimates are not always sharp. Furthermore, these experiments will show these terms 
to be computationally impractical. 

Moving now to the flux formulation [equation (80)], define the vectors u = (u(xo,t), ■■■ ,u(x n ,t)] T , and (j) = 
[u x (xo,t), ■ ■ ■ , u x (x n , t)] T , to be the projected values of the exact solution and derivative, defined in the left domain 
at gridpoints x. Similarly, define the vectors v, and in the right domain. With these definitions, the continuous 
problem (30) can be represented in vector nomenclature as 


u t + au x - €(j) x 

— Loo\ 
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- Vi] 

e i- 
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K 1 e i+ 
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Rii (ip - 

<t>]i 

Vf 1 e i+ 


with suitable initial and boundary data, and the source term F(x, f) = 0. Note that the interface terms [ui — vf\ and 
[4> — if\i are identically zero for the exact solution, but have been added to the right hand side of equation (126) for 
added clarity. 

Substituting the discrete accuracy relations u x = 2 ?/u + T e ii , v x = 2 ?/v + T e i r , f> x = Dif) + T^, and 
Vtc = 'D r f> + T .0 into equation (126), subtracting the result from V x equation (80), and defining the semi-discrete 
errors Hi = U — u, E r = V — v, and ©i = (f> — tp, © r = tj> — ip, in the left and right domains, respectively, 
produces an evolution equation for the error. Specifically, the collocation error of the coupled problem (168) written 
in flux form is 
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(127) 
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K 1 e i+ 
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Vp 1 e i+ 



eT e i r 


S P (a:,0) = 0. 

Equation (127) is not in the form of a semi-discrete PDE. Thus, the accuracy of the flux form solution can not be 
inferred directly from Theorem (11). To assess the accuracy, equation (127) must be transformed from flux to primal 
form [For details on how to relate the two forms, see proof of Theorem (8)]. 

Solving equation (127) for ©i and © r yields the expressions 

+ h'Pl 1 {-^lo[(S0i - (Sj)J e i- + An[(0;) i - (0 r )J ei } 

+ 7'P r r 1 {7?io[(S;)j - (S;)J e i+ + i?n[(0 r ) i - (0;)J e i+ } 
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Rewriting equation (127) in terms of the difference (0;) i — (0^ yields 


(©/), - (0r)i = <*[(3^ - (S r )J + c 2 [(ASOli - (V r S r )|i] - C 2 [(T elI )|i - (T elr )|i] (129) 

with the constants ci and c 2 defined in equation (101). Substituting equations (128) and (129) back into (127) yields 
the error evolution equation of the LDG flux form in terms of primal form variables. 
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Four distinct types of discretization error appear in equation (130). The first derivative truncation error vector: T e i 
matches that found in equation (124). The second derivative truncation vectors (DT e i + T^,) are dissimilar to those 
found in (124), but will have similar leading order truncation terms. The final two truncation vectors result from the 
imposition of penalties and have the same form as those found in the primal form (124). Thus, based on this analysis as 
well as the truncation analysis performed on the primal formulation, both should have the same asymptotic behavior 
with respect to grid refinement. 


6 Identification of schemes 

Specific values are assigned for the eight parameters used in the primal [equation (46)] and flux [equation (80)] 
formulations. Three popular schemes are identified, including two DGFEM and one strong form HOFDM. Each is 
shown to satisfy their respective stability constraints. 

6.1 Carpenter, Nordstrom, Gottlieb [13] 

We begin with the scheme proposed by Carpenter et. al (CNG) [13] which is a special case of the primal formulation 
presented in equation (46). This scheme was originally derived in strong form in the context of HOFDM. This scheme 
[See reference [13], equation (7).] can be reproduced by assigning the following specific values to the penalty parameters; 

ho = hi = Do = rn =0 ; loo = cji ; Zoi = ecr 2 ; roo = 03 ; Hu = £ 0-4 ; (131) 

Substituting (131) into equations (47, • • -, 50) yields the conditions 

01 = 03 + a ; [<t 2 = 04 — 1 ] e ( 132 ) 

0 < [(a + (3)e\ ; ]e(-a + /3 )] 2 < [(a + /3)e ] 2 ; 01 < - - + 7 ^ 7 ] (133) 

Further, the stability conditions obtained in equation (132) and (133) are identical to those found in the original CNG 
work [e.g. equation (8) in reference [13]]. 

6.2 Baumann, Oden [7] 

The scheme proposed by Baumann and Oden [7] is a special case of the primal formulation presented in equation 
(46). This scheme was originally developed in the weak form in the context of DGFEM. Shu presents in reference 
[35], an illustrative comparison between this scheme and the LDG scheme for the strictly parabolic case. Therein the 
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Baumann- Oden scheme is presented in equations (5.1) and (5.2) and can be reproduced from equation (46) with the 
■values 

Zoo = r oo = hi = rn = 0 ; Zoi = — Zio = — ^oi = + ri o = + (134) 

Taking the liberty to extend Shu’s definition to include the possibility of convection terms, we define the Baumann- 
Oden scheme by assigning the following values to the coefficients in equation (46): 

Zoo = r 00 + a; hi = m = 0 ; l 01 = +(/3 - |); ho = + |); r 0 i = +(/3 + |); ri 0 = -(/3 - |) (135) 

Substituting (135) into equations (47, • • •, 50) yields the conditions 

0 < [(a + /?)e] ; [e(-a + /3)] 2 < [(a + fi)e} 2 ; Z 00 < | (136) 

6.3 Local Discontinuous Galerkin [14] 

The scheme proposed by Cockburn and Shu [14] is special case of the flux formulation presented in equation (80). 
Once again, Shu presents in reference [35], the LDG scheme written strictly for the parabolic case. Therein the LDG 
scheme is presented in equations (4.1), (4.2), and (4.3) and can be reproduced from equation 80 with the values 

Loo = Roo = Ln = Rn = 0 ; R 01 = — — ; A 10 = — R 01 = L Rio = + ^ (137) 

Taking the liberty to extend Shu’s definition to include the possibility of convection terms, we define the collocation 
LDG scheme by assigning the following values to the coefficients in equation (80): 

Loo = -Roo + a; Ln = Rn = 0 ; A01 = +(/? — ^io = — + 9)’ R01 = +(/? + 9); Rio = —{P — ( 138 ) 

Substituting (138) into equations (47, • • •, 50) yields the conditions 

0 < [{a + f5)e] ; [e(— a + (3)\ < [(a + /3)e] ; Roo < ^ (139) 

7 Numerical experiments 

7.1 Test Problems 

7.1.1 The One-Way Wave Equation 

The linear one-way wave equation is used to study the stability and accuracy of the new formulations in the absence 
of diffusion terms. The functional form is 

U t + aU x = 0 |z| < 1 ,0 < t < 0.1 (140) 

with the exact solution 

U(x,t) = cos[— 2 7T (x — t )] (141) 

The initial and boundary data coincide with the exact solution for all time. 

The wave equation is used to establish a baseline accuracy for each method, before moving on to the parabolic 
equation. This step is essential to establish the effects of the newly added terms that penalize the diffusive terms. 
Specifically, the new additional terms could destroy the accuracy of the advection term, not just diffusion. Comparing 
the two results allows us to determine the source of error. 

This study of the first order wave equation raises a subtle point regarding imposition of discrete interface conditions; 
specifically, that derivative fluxes can be penalized between elements, despite the fact that wellposedness requires only 
one interface condition. Thus, the diffusive interface penalty can be thought of as an alternative discrete stencil, 
one that enforces stronger first derivative continuity across the interface. An a priori stability estimate and accurate 
interface data guarantee the validity of these penalties. 
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Figure 4: Schematic depicting two periods of the recursive element pattern used in the nonuniform grid accuracy 
studies. 

7.1.2 Linear Burgers’ Equations 

The linear Burgers’ equation tests a combination of the advection and diffusion terms. Five distinct values of the 
diffusion parameter e are used to deduce the effects of the penalty as the problem becomes diffusion dominated. The 
functional form is 

Ut. T aU x = eU xx \x\ < 1 , 0 < t < 0.1 (142) 

with the exact solution 

U(x,t) = exp[— eb 2 t] cos[— 2 tt (x — t)] (143) 

with initial and boundary data that coincides with that of the exact solution. The equation parameters used are a = 2 
and b = 2 7r in all cases. Five values of the parameter e are used: e = 1 , |, 0.1 and 0.01. 


7.2 Test Schemes 

Two finite difference and four spectral operators are tested using the linear wave and linear Burgers’ equation. The 
finite difference operators are the 4 th and 6 th order central difference operators using SBP closures at the boundaries. 
The boundary closure formulas are those reported in reference [10]. The spectral operators are conventional collocation 
operators defined using the Gauss-Lobatto-Chebyshev nodal points. Polynomial orders in the range 2 < p < 5 are 
used (corresponding to elements having three to six nodal points). The penalties are implemented such that the scheme 
is equivalent to a Legendre collocation scheme (see Carpenter and Gottlieb [12]). The time advancement scheme used 
in all cases is the five-stage fourth-order Runge-Kutta (RK) scheme [9] . The timestep is chosen such that the temporal 
error estimate is independent of further timestep reduction. The magnitude of the temporal error is approximately 
10“ , and is obtained by comparing the fourth-order solution with an embedded third-order solution. Thus, the 
spatial error is the dominant component of error in the simulations. 


7.3 Test Matrix 

An extensive test matrix is used to distinguish different properties of the new formulations. The matrix includes 1) four 
geometry definitions, and 2) multiple values of three distinct interface penalty parameters. The geometric definitions 
included all permutations of uniform and nonunifornr element distributions on both periodic and nonperiodic intervals, 
specifically; 1) uniform-periodic, 2) nonuniform-periodic, 3) uniform-nonperiodic, and 4) nonuniform-nonperiodic. The 
nonunifornr grids are generated using elements that alternate in size with a ratio of 2 : 1. Figure 4 shows a schematic 
depicting two periods of the 2 : 1 nonunifornr grid used in the study. 

An a priori study was performed, comparing the accuracy of the methods on each of the four different geometries. 
The results from one of these studies are included in Appendix (E). The most challenging test cases (hyperbolic 
or parabolic) were those that were run on the nonuniform grid, with nonperiodic boundary conditions. Thus, the 
nonunifornr grid with nonperiodic boundary conditions is used almost exclusively in this work. 

The interface penalties are parametrized using three independent variables as (here written for the primal formu- 
lation), 

^oo = ^(1 - «); r oo = |(-1 - a); a > 0 (144) 
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(145) 


hi — +{(3 — -); ho — — {(3 + -); r m — +(/3 + -); rqo — — ((3 — — ) 

hi = rn = 7; 7 > 0 (146) 

The parameter a adjusts the contribution of the Dirichlet penalty between the left and right states. The special value 
a = 0 admits no dissipation at the interface, thus producing a skew-symmetric matrix in the "P-nornr. The value 
a = 1 produces a fully upwind flux. Two values of the parameter a are tested: a = 0 and a = 1. 

The parameters (3 and 7 adjust the penalty on the derivative fluxes. For (3 the relationship between the values of 
Z01 ) ^10) r 0 i, and 7*10, greatly simplifies the energy estimate by producing the conditions t\ = <2 = 0 in the estimate. 
This combination of parameters produces no contribution to the energy estimate. [See equation (62).] Note that the 
value of (3 can be chosen independent from e, and can be nonzero for values of e = 0. Two values of the parameter 
(3 are tested: (3 = 0 and (3 = j. The parameter 7 adjusts the contribution from the derivative fluxes and produces 
contributions to the energy estimate that are strictly dissipative. Four values of 7 are tested: 7 = 0, 7 = 10 -4 , 
7 = 1CT 3 and 7 = 1CT 2 . 

All studies include a comparison of the convergence rate of each scheme. A grid refinement exercise involving 
approximately 40 grid densities is performed for each set of parameters. The convergence rate is determined by fitting 
a “least-squares” curve through the error data. 

7.4 Results 

7.4.1 Finite Difference: Wave and Burgers’ Equation 

Figure 5 compares the convergence behavior of the 4 th - and 6 t/l -order finite difference schemes. The model problem 
used is the linear one-way wave equation. Shown are the results obtained on a nonuniform grid with nonperiodic 
boundary conditions. Block- norm boundary closures are used that are 3 rd - and 5 t/l -order accurate [10], respectively. 
Seven sets of interface parameters are compared, and in all cases, design order convergence is observed (4 th - and 
6 t/l -order, respectively). Note that imposing the penalty on the derivative-transpose terms (in, and rn) does not 
decrease the formal accuracy of the schemes. 

Figure 6 compares the convergence behavior of the the A th and 6 t/l -order schemes on the linear Burgers’ equation. 
Both the flux (LDG) and the primal (B-O) forms of the penalty are presented in this study. Note that the parameter 
(3 is not scaled with the parameter e. In all cases, design order is achieved, implying that increased levels of interface 
coupling has little impact on formal solution accuracy. Two other values of the physical viscosity e were also simulated 
(but not shown) to assess the impact of the interface penalties. Similar trends were observed at other values of e. 

An anomalous behavior is observed in the convergence of the 4 tft- -order simulations. While the slope remains close 
to 4, a deflection in convergence behavior is observed at an error of approximately 10 -8 . The exact cause of this 
behavior is not known. 

All simulations were performed with an explicit RK scheme, and the efficiency of the runs is strongly influenced by 
spurious eigenvalues generated by the interface coupling terms. For the parameters tested, little sensitivity is observed 
for the a parameter. A significant decrease in the max-CFL is observed when large values of the (3 parameter are used. 
Even small values of the 7 parameter rendered the explicit method almost useless because of the magnitude of spurious 
eigenvalues produced. An implicit method would most certainly be needed if the 7 parameter is used extensively to 
stabilize the interfaces! 

Finally, both uniform and nonuniform grids were used in this study, as were periodic and nonperiodic geometries. 
Neither permutation changed in any way the general convergence behavior of the schemes presented in figures 5-6. 

7.4.2 Spectral Collocation: One-Way Wave Equation 

The same a , (3 , 7 parameter studies performed for the finite difference cases are repeated using the spectral collocation 
methods. Polynomial orders ranging from p = 2 to p = 5 are compared on a nonuniform grid with nonperiodic 
boundaries. In addition to graphical comparisons, a numerically determined convergence rate is used for direct tabular 
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comparisons. Three different measures of error are used in the tabular comparisons: 1) the L 2-norm of the solution 
error, 2) the Loo-norm, and 3) the integral norm, formed using the integration matrix V . 

Table 2, shows the convergence results for four permutations of a and /3, and polynomial orders p = 2 - p = 5, 
at values of 7 in the range [0.0 < 7 < 0.025]. Three distinct convergence rates emerge in this study. We begin by 
comparing the cases run with 7 = 0. The case a = 1, j3 = 0 converges at the optimal rate of p+ 1 for all polynomial 
orders. Changing the value to (3 = j reduces the formal convergence rate to p. Next, imposing a nonzero value of the 
parameter 7 decreases the convergence rate of the optimal case ( a = l,/3 = 0) to what appears to b e p+|. Strangely 
enough, imposing nonzero values for 7 increases the accuracy of the suboptimal cases for which /3 = j. 

Once again, all these studies were run with an explicit RK scheme that is subject to a finite maximum stability 
condition. The interface penalties have a profound influence on the placement of the maximum eigenvalues, and 
therefore the maximum allowable timestep [(A t) max ] . The parameter (3 decreased the (A t) max significantly for 0( 1) 
values. The parameter 7 rendered the explicit time advancement scheme almost useless, by reducing the (A t) max by 
several orders of magnitude. 

7.4.3 Spectral Collocation: Linear Burgers’ Equation 

Two studies are performed on the linear Burgers’ equation: one using the primal formulation (e.g. Baumann-Oden), 
and one using the flux formulation (e.g. LDG). We begin with a study comparing the influence of the equation 
parameter e on the convergence rate of the Baumann-Oden scheme. 

Table 3 compares the convergence of the B-0 scheme for values of the physical parameter 
e = 1.00 , 0.10 , 0.01. (The I/2-norm of the derivative error is included in the results along with the Z/2-norm. the 
Loo-norm, and the V integral norm.) The cases with diffusive parameter e = 1.0, are dominated by the parabolic 
nature of the equations, while as e becomes smaller, the character of the equation becomes increasingly “hyperbolic”. 

Suboptimal convergence is generally observed for all polynomial orders and schemes. The exceptional cases are 
when the Dirichlet boundary and interface terms are fully upwinded (a = 1) and the Neumann derivative terms are 
averaged (/3 = 0). For this set of parameters, polynomials of odd order (p =: 3,5) converge at an optimal rate of 
p + 1 for the solution and p for the derivative. This result is consistent with that shown by Shu [35] for the strictly 
parabolic case using DGFEM. As e decreases from the baseline e = 1, the convergence rate increases (unexpectedly), 
for all schemes and all polynomial orders with the exception of those that are already optimal. 

The convergence behavior of the B-0 formulation is unpredictable. In contrast, when the linear Burgers’ equation 
is simulated with the LDG scheme, optimal convergence is achieved for every combination of physical/numerical 
parameters. Table 4 shows the convergence behavior of the LDG scheme for the same four parameter settings, and 
four polynomial orders. The LDG scheme converges for the solution at a rate of p+ 1, while the derivative converges at 
a rate p for all polynomial orders p. Little influence is observed when the values of the parameters a or (3 are adjusted. 
These results for the LDG scheme are in slight disagreement with those presented by Shu [35] for the DGFEM and a 
strictly parabolic equation. There it is shown that optimal convergence is only achieved for the specific values of the 
parameter (3 = ± 

The formal accuracy of the LDG scheme is clearly superior to that of the B-0 formulation. Figures 7, 8, 9, and 
10 compare the absolute accuracy of both formulations, all permutations of a and /3, and the physical viscosity e = 1. 
The LDG formulation significantly outperforms the B-0 scheme for even order polynomials, but is less advantageous 
for odd order polynomials (for which the B-0 scheme is optimally convergent). Both the B-0 and LDG classes of 
schemes, exhibit little sensitivity to the parameters a = 0, 1 and (3 = 0/4, 1/4. This result is surprising for the 
B-0 schemes, considering that some are suboptimal for various parameter combinations. Note, however, that absolute 
accuracy depends on both the order of and size of the truncation terms. Many of the lower-order combinations have 
smaller truncation terms, which results in better accuracy at coarse tolerances. 

Figure 11 combines all the polynomial orders in to one plot. For clarity, only the cases with interface parameters 
a = 1 and f3 = 0, are presented. Recall that this parameter combination results in an upwind Dirichlet flux and a 
centered Neumann flux at the interfaces, and resulted in optimal convergence for odd-orders with the B-0 scheme. 
The LDG scheme is significantly more accurate than the B-0 scheme for even order polynomials, while they are nearly 
indistinguishable for odd order polynomials. Note that the convergence behavior of the LDG formulation is dual valued 
for even order polynomials. Specifically, the simulations with an even number of elements superconverges initially, 
while those with an odd number of elements converge as expected. 
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8 Discussion 


The present work unifies existing FEM methods into a common strong form framework, and then generalizes that 
framework by including additional interface penalty terms. The result is a “plethora” of adjustable parameters that 
generate schemes with remarkably diverse numerical properties. We have provided some broad generalities regarding 
the choice of parameters, e.g. the accuracy virtues of the LDG flux formulation. Many readers will want concise guid- 
ance regarding optimal parameter choices. Given the subjective nature of optimality, however, we have intentionally 
avoided advocating specific parameter values. The following examples illustrate the reasons for our choice. First, a 
practitioner using an implicit temporal algorithm could utilize all first derivative terms, as well as second derivative 
penalty terms to enrich interface continuity. This is in contrast to the general conclusions of the present study, i.e. 
that second derivative terms are never advisable if explicit temporal algorithms are used. Second, specific choices 
of penalty parameters in the flux formulation result in nearest neighbor stencil compactness, a property that some 
(but not all) practitioners find desirable. Third, nonlinearity is the ultimate test of scheme accuracy and robustness. 
Excluding portions of the parameter space based on the present studies in the absence of limiters is premature. 

The present derivations and analysis, focus entirely on the scalar, constant coefficient, advection-diffusion equation. 
This lack of generality is an obvious shortcoming for most practitioners, given their desire to solve systems of nonlinear 
equations. Straightforward generalization to systems of equations is possible for the penalty terms that originate in 
the multi-domain work of Carpenter, Nordstrom, Gottlieb[13]. (See the work of Nordstrom, and Carpenter [30] for 
details.) Generalization (FEM reverse engineering) of the penalty terms involving V T to systems of equations is the 
focus of ongoing work. 


9 Conclusions 

Summation-by-parts (SBP) operators are discrete derivative operators that mimic the integration-by-parts property 
of the continuous operator. A diverse set of discretization operators can be cast in the SBP convention, including 
liigh-order finite difference, (central, upwind, compact, DRP) and spectral collocation operators. Summation-by-parts 
operators share a close resemblance with finite element methods. In fact, for certain polynomial methods, the (strong 
form: collocation) SBP operators can be shown to be equivalent to (weak form: integral) FEM methods. The great 
similarity between the SBP and FEM operators [especially Discontinuous Galerkin FEM (DGFEM)], leads to a fertile 
exchange of ideas between the two approaches. 

Two new collocation approaches are developed to couple multiple SBP domains. The two approaches are essentially 
the strong form equivalent of the primal form method developed by Baumann and Oden, [7] , and the flux form LDG 
methods developed by Cockburn and Shu [14]. Both techniques are shown to be L 2 and pointwise stable, and 
conservative, while (if applied carefully) maintaining the original design accuracy of the adjoining domains. The new 
approaches generalize the original multi-domain penalty technique of Carpenter et al. [13] . 

The Baumann and Oden collocation approach is extended to allow construction of interface penalties involving 
second derivatives. Several free parameters are available in the new approach, that can be adjusted to ensure L 2 and 
pointwise stability of the interface and two adjoining subdomains. The general structure of the new approach enables 
(in principle) penalties of derivatives higher than second to be constructed. Similar ideas can be used to construct 
generalized penalties based on the LDG collocation approach. 

Most penalty combinations are theoretically advantageous because they enhance the interface “continuity” . There 
are, however potential costs. Penalizing data that is of inadequate accuracy can degrade the formal global accuracy of 
the method. Penalties built from second and higher derivatives are ill advised for this reason. A further potential cost is 
the degree to which some penalties increase the condition number of the resulting discretization matrix, making them 
unacceptably expensive for explicit time advancement. Extensive numerical experiments on the advection and the 
advection diffusion (linear Burgers’) equation verify the theoretical predictions, although in some cases the theoretical 
predictions are not sharp. 
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A FEM Stability 


Although arbitrary interface closures are admissible for equations (9) or (10), not all choices result in stable or accurate 
formulations. An interface parametrization and derivation of the L 2 stability of the LDG method is now presented. 
As we shall see (section III), this stability analysis closely resembles that used to define admissible strong form 
discretizations. 

The LDG energy is obtained from equation (3) by choosing the test functions to be v = u and w = q, resulting in 

In u t u dx = - f n qu x dx + Er [[<H] 

Jn q qdx = -f Q uq x dx + Er [Ml 

Adding the two “energies” yields after manipulation the equation 

IJi.In u2dx + fn ^ 2 dx = Er [[«« +uq~ <?«]] 

To this point, the functional form of the interface variables u and q, remains undefined. Their choice, however, greatly 
influences the properties of the resultant scheme. In general, one would like to choose these fluxes such that the global 
scheme is stable, conservative and exhibits optimal accuracy. 

The fluxes u and q are generally defined to be conservative if they are single- valued on all interfaces T. Conservative 
fluxes that telescope across domain interfaces, eliminate many interface terms in proofs of stability and accuracy. For 
example, inserting v = 1 and w = 1 into equation (6) produces the trivial condition f n ut.dx = f n q = 0, for 
conservative fluxes. 

Interface accuracy can only be maintained if the interface fluxes are constructed from accurate interface data. 
Since all solutions are double valued on the interfaces, natural candidates are averages and differences of the interface 
variables. This in mind, one might choose a conservative interface flux of the form 

u = {u} + Gi[[u]] + C 2 [[q}\ n . q , 

? = {<?} + C 3 [[u}} + C 4 [[q ]] 

Substituting these fluxes into the energy estimate given in equation (148), using the relation [[£77]] = {C} [[??]] + [[C]] {dA 
and the conservation conditions [[it]] = [[<)]] = 0 produces 

~[ u 2 dx + [ q 2 dx = Y, (Cl [[«]] + C 2 [[</]]) [[«]] + (C 3 [[«]] + C 4 [[g]]) [[q]} (150) 


(147) 

(148) 


The (RHS) equation (150) is negative (implying stability) if all the conditions 

Ci <0, C 4 < 0, [C-! C 4 - f(C 2 + C 3 f] > 0 (151) 

are satisfied. The precise values of Ci , • • • , C 4 are not important for this discussion, just the functional form of the 
fluxes. 

The primal formulation given by equation (9) produces a slightly different energy estimate, because stability is 
proved in different variables. Choosing the test function to be v = u in equation (9) produces the energy equation 

\liiSn u2dx + fn u l dx = +' du x ~ uu x }\ . (152) 


Once again, one might choose a conservative interface flux of the form 


u = {u} + ci [[uj] 

u x = {u x } + C3 [[«]] 

Substituting these primal fluxes into equation (152) produces 


C2 [[Az]] 
C4 [[«*]] 


1 d_ 

2 dt 


u 2 dx + 



E ( £ 1 [Ml + C 2 [[«*]]) [[«]] + (C3 [M] + C4 [[«*]]) [[«*]] 


(153) 


(154) 
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The methods are stable if the stability conditions 


are satisfied. 


ci <0, C4 < 0, [ci C4 


\{ c 2 + c 3 ) 2 ] > 0 


(155) 


Remark. The stability conditions of the two energy estimates are identical with the exception that different diffusive 
variables appear in the proofs. 

Remark. The admissible parameters for the interface fluxes u and q (or u x ) cannot be chosen independently. 
Remark. A reasonable (and somewhat more clever) person might choose a conservative interface flux of the form 


u = {u} + ci [[«]] + c 2o [[q]} + c 2 b [[«x]] 

u x = {u x } + c 3 [[«]] + c 4a [[q]} + C4,b [[u x ]]. 


(156) 


Numerous methods, similar in type to the LDG and IPFEM methods have recently been developed. They differ 
primarily in the functional form of the interface fluxes. Arnold et al. [2, 3] presented a framework that unifies nearly 
all the DG schemes that appear in the literature. A subset of the many schemes presented in Arnold et al. [3] is 
included in Table 1. 


Method 

u 

q 

Bassi-Rebay [5] 

M 

{q} 

Brezzi et al [8] 

M 

{q}-a r [[u]\ 

LDG [14] 

M - PIM 

{9} - otj[[u]\ + P[[u]} 

IP [16] 

M 

{u x } - CXj [[«]] 

Bassi et al. [6] 

M 

{«x} - a r [[u]} 

Baumann- Oden [7] 

M + n ■ [[a]] 

{u x } 

NIPG [34] 

M + n ■ [[u]] 

{u x } - aj[[u]\ 


Table 1: Numerical fluxes for various DG schemes. 


Note that the methods fall into two distinct categories based on how they form the interface flux q; the first three 
penalize the flux variable q , while the last four penalize the derivative of the solution u x . 


B The continuous problem 

In this paper we will consider interface conditions between subdomains. However, interface conditions are closely 
related to boundary conditions; therefore, we start with the single domain problem. 


B.l The continuous single domain problem 


We begin with a discussion of the wellposedness of the linear advection-diffusion problem (Burgers’ equation), and 
present a simple proof. 

For brevity, we reintroduce the linear Burgers’ equation [originally presented in equation (30)], as 


c t + a c x 

= £C XX + F(x,t) 

O 

Al 

•SO 

,-l < X < 1, 

c(x, 0) 

= f{x) 

O 

II 

-to 

,-1 < X < 1, 

7c(— l,f) - ec x (—l, t) 

= 9-i(t) = L -1 (c) 

0 

Al 

-to 

,x = -1, 

C c ( + 1’ 1) + ec x( + li t) 

= g+i{t) = l + 1 {c) 

0 

Al 

-to 

, X = +1, 


(157) 


with constants a and e (0 < e). Arbitrary values of the parameters 7 and ( cannot be used, but they can be determined 
from the proof of wellposedness. 


Theorem 1 Burgers’ equation (15 7) is strongly wellposed with boundary conditions satisfying the conditions 


0 < a + 2 ( ; 0 < — a -t~ 2 7 
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Proof : An energy estimate for equation equation (157) is obtained by multiplying by c and integrating over the 
domain — 1 < x < 1 and takes the form 

|l|c||n = -2e|| c x \\l + 2\\cF\\ n + [-ac 2 + (158) 


Solving the boundary conditions in equation (30) for ec x yields the relations 


ec x |_i = (— ff-i + 7 c -i) ; 

ec x \ +1 = (+3+i-Cc+i) 

(159) 

which upon substitution into equation (158) yields 



ilMln = - 2e|| c x f n + 2\\cF\\ n 

+ 2 c+i 3 +i — ( 2£ + a)c+i 2 

+ 2 c_i 3 _i - ( 2y - a)c_i 2 

(160) 

Expressions of the form 



2AB = -(y/lA-^B) 

+ fA 2 + jB 2 ; 0 < £ 

(161) 


can be used to reduce equation (160) to the form 


mMn + 2 e \\ c x\\n = (Ml - ( 2 C + a + £+i)c+i 2 - (2y - a + £-i)c-i 2 

+ fll^Hn + 777 5+l2 + (162) 

7? F Hn “ - (-v/^Tc-i - -jL=g-!) 

Grouping like terms results in 

i{e- ?t ||c||^} + e - ^ {2 e ||c x ||q + BC} = e~^ [Data - PosDef } (163) 

with 

BC = + (2£ + a + £ + i)c+i 2 + (2q — a + £_i)c_i 2 

Data = + ±\\F\\ 2 n + ^g +1 2 + ( 164 ) 

PosDef = + || Vic - ^-F|| 2 + (y^C+i - -±=g +1 f + (y^C-i - -^= 5 _!) 2 

The wellposedness conditions for the boundary data, are determined by the requirement that the BC term be positive 
definite. The fact that £_i and £+i are arbitrary positive constants, imposes the following constraints on the boundary 
data 

0 < a + 2 ([ ; 0 < —a + 2 7 . (165) 


The Data term is problem specific, while the PosDef term is positive definite by construction. Integrating equation 
(164) in time and imposing initial conditions yields 

\\c(x,T)\\ 2 n + / 0 T e-«(‘- T ) {2 e ||c x ||q + BC}dt. = e+« T {||/(:r)||^ + e"«‘ [Data - PosDef]dt } (166) 

from which a strongly stable estimate of the form 

\\4 2 a + r ||c || 2 v dt < K c e^{\\f\\l + A||F|| 2 + \\g\\ 2 )dt}, (167) 

Jo Jo 


immediately follows. □ 

Remark. A sharper proof is given in great detail in the work of Hesthaven and Gottlieb [20] and will not be repeated 
here. 
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B.2 The continuous multiple domain problem 


coupled problems are 


it the domain [— 1, 1] into 

[—1,0] and [0, 1] and focus 

on the interface problem , 

u t + au x 

= eu xx + Fi(x,t) 

,t> 0 

, — 1 < x < 0, 

u 

= fi(x ) 

O 

II 

-to 

,-l < x < 0, 

7u(-l,f) - eu x (-l,t) 

= 9-i{t) = L-i(u) 

,t> 0 

,x = -1, 

1? 

JS, 

c 

= 0 

,t> 0 

,x = 0, 

V t + av x 

= ev xx + F r (x,t) 

,t> 0 

,0 < x < +1, 

V 

= fr(x) 

o 

II 

-to 

,0 < x < +1, 

o 

i 

= 0 

,t> 0 

,x = 0, 

C ^(+1; t) + e Ue(+1, t) 

+ 

II 

+ 

II 

O 

Al 

-to 

, X = +1, 


( 168 ) 


respectively. The variables in the left [—1, 0] and right [0, +1] domains are u and v, respectively. The coupling between 
the two domains in equation (168) is given by the operator L(x = 0). 


B.2.1 Wellposedness 


Subtracting (30) from (168) and assuming linearity of L o, produces the equation 


{5u) t + a(6u) x 

= < 5u )xx 


,t> 0 ,-l < x < 0, 

{6u) 

= 0 


, t = 0 , — 1 < x < 0, 

7 (<^u)(-l, t) - e(6u) x (- 1, t) 

= 0 = 

L -l(( su )) 

,t> 0 , x = —1, 

L 0 ((Su) - (Sv)) 

= 0 


IV 

o 

H 

II 

O 

(5v) t + a(5v) x 

= < Sv )xx 


,t> 0 , 0 < a; < +1, 

(Sv) 

= 0 


,t = 0 , 0 < x < +1, 

L 0 ((Sv) - ( Su )) 

= 0 


IV 

o 

H 

II 

O 

C (<^)(+M) + e (<^)x(+M) 

= 0 = 

L +i((5v)) 

,t> 0 , x = +1, 


(169) 


where ( Su ) = u — c and (5v) = v — c describe the evolution of the differences between the one and two domain 
solutions. The solutions to equation (169) depends on their boundary conditions imposed at the interface point x = 0. 
It must be shown that imposition of this boundary condition does not affect the wellposedness. 


Theorem 2 If Theorem 1 is used to provide exterior boundary conditions, and the interface conditions 

f Ui-Vi \ _ f ( Su)i - ( Sv)i \ _ „ 

V ) V [(&*)*- OH*] I j ) 


are used, then the semi-discrete approximation on two domains (168) is strongly wellposed. 


Proof : The energy method applied to equation (169) leads to 

JHlIWIln + \\(Sv)f n } = -2 C ||(5 U )X + [(6u)a(6u)-2e(6u)(6u) x r x z°- 1 

- 2 e ll(<HJIn + \{Sv)a(5v) - 2e(5v)(5v) x \ x Zo 


(170) 


The stability of the single domain problem implies that the boundary terms at x = ±1 are negative semidefinite 
with the boundary operators given by Theorem 1. The terms arising at the interface x = 0 that contribute to the 
estimate are 


[(Su)a(Su) - 2e(Su)(5u) x ] x 0 + [(<5v)a(<5v) - 2e(5v)(5v) x ] x = 0 



/ (Su)i + ( Sv)i \ 

T 

( 0 

a 0 

— e 

\ 

/ (5u)i + ( Sv)i \ 

1 

(Su)i - (5v)i 


a 

0 -e 

0 


(Su)i - (Sv)i 

2 

(( Su) x + (<5v)x)|i 


0 

-e 0 

0 


N* + 


\ ((Su) x - (Sv) x )ji ) 



0 0 

0 

/ 

V [(<M* - (H*]li / 


( 171 ) 
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The eigenvalues of this symmetric matrix are (±a± \/a 2 + 4e 2 )/2 which includes both positive and negative values. A 
negative semi-definite contribution for the interface terms can only be achieved by choosing the boundary conditions 
(Su)i — (Sv)i = 0 and [(<5«) x — (<5v)z]lj = 0, the desired result. □ 

Remark. The two domain problem (168) is strongly wellposed in the sense that the solutions can be estimated in 
terms of the data in the corresponding one domain problem (30). 


C Positive Semi-Definite Norms 

A small portion of the diffusion term ||2?U||.p arising at the boundary of a domain (element) can be moved to the 
(RHS) to help in the proof of stability. This “borrowing” must be accomplished without destroying the definiteness 
of the remaining term ||2?U||p where 

V = V - ae i+ e i+ T (172) 

with ei + = [1, 0, ..., 0, 0]. What remains to be shown is the range of a for which the new matrix V is positive 
semi-definite. 


Theorem 3 Assume that the original matrix V is symmetric positive definite. A positive constant a can he subtracted 
from the Vn position, without destroying the definiteness, if the value of the parameter a is in the range 


0 < a < 


l 

e i+ T V- 1 


(173) 


Proof : The new matrix P is still symmetric, and the definiteness of a symmetric matrices is completely determined 
by its eigenvalues. Thus, if the eigenvalues A j of V satisfy 0 < A,, then the matrix is positive semi-definite (positive 
definite for the strict inequality). The eigenvalues A j are identical to the roots of the characteristic polynomial. 
Accounting for the special form of the matrix V, the characteristic polynomial is: 


detfiP - XI) 


Expanding the determinant yields 


Pn - ol - A P 12 
Pv2 P22 — A 


Pin 

P2N 


Pin P2N • • • Pnn — A 


(-1) JV (A) JV + {—l) N ~ 1 Trace\P]{X) N ~ 1 + ••• + detfV) = 0. 


(174) 


(175) 


A zero root for equation (175) can be obtained by requiring the term det(V) = 0, a condition that allows a A to be 

factored out of the expression. Expanding the expression det{ V) = 0 produces the following condition 

det{V) = det(V) — adet(Mn) = 0 (176) 

where Mu is the (n — 1) x (n — 1) minor of the original matrix P obtained by simultaneously eliminating row and 
column one. Solving equation (176) for a (and using Cramer’s rule for the element Vfi) yields the desired result. 

— det(V) — 1 — 1 1177) 

a — det(Mu) ~ — el T V 1 e,_ \ L ‘ ‘ ) 


a. 

Remark. The result of Theorem 3 is trivial in the case of the diagonal norm. 

Remark. P The corresponding proof that accounts for borrowing a from both the Vn and Vnn terms has not 
been constructed. 
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D Pointwise Boundedness 


We begin with the Sobolev inequality by presenting without proof the following continuous and discrete theorems. 
Theorem 4 Let u € C 1 (0 < x < /). For every e > 0 

NIL < r 1 + + e||«x || 2 (178) 

Proof : See Gustafsson, Kreiss, and Oliger, [19], Lemma 9.3.1, page 377-378, for the proof. □ 

Theorem 5 For any grid function U on (0 < x < l ), and every e > 0, we have 

HUIlL < (e- 1 + r^liuf + ell^uf (179) 

with the discrete derivative operator T>\ being defined as 

-1 1 0 
-1 1 0 


-1 1 0 
-1 1 

Proof : See Gustafsson, Kreiss, and Oliger, [19], Lemma 11.2.1, page 459-460, for the proof. □ 

Remark. Note that the derivative matrix T>i is of dimension (n-lx n ), that each row is linearly independent, and 
the row-sums are zero. The rank of the matrix is therefore n — 1, and the null space is the constant vector. 

The energy method applied to the linear Burgers’ equation [see equation (30)], provides ’’integral” bounds on 
the solution, and the discrete first derivative. Specifically, the quantities ||U|jp and ||PU||p are bounded for finite 
times. The original discrete Sobolev inequality is proved using the first order derivative operator T>\. Thus, to use 
this estimate for pointwise boundedness, we first must establish that P-boundedness of ||U||^ and ||2?U||^,, implies 
P-boundedness of ||U||^ and ||2?iU||j. 

The norms V, and X are of full rank n. Thus, that ||U||p implies HU^ is immediate. Because of the inherent 
difficulties in dealing with rank deficient matrices, the derivative condition is more subtle. 




Theorem 6 For any grid function U on (0 < Xj < l) ; j = 1, n, and a consistent derivative operator V of rank n — 1 
(the constant vector is the only null vector), then V -boundedness of ||PU||p, implies V -boundedness of ||2?iU||j. 

Proof. Any consistent , first derivative operator, consisting of n distinct stencils for the n grid points on the interval 
Xj, can be represented as linear combinations of the n — 1 rows of the matrix T>\. To illustrate this, imagine row i of 
T> written as 


P | i — [b> * * * i b? di,i— , • • • , d,i /l+ k r , 0, * * • , 0] — [0, • * * , 0, Tflij—ki , * * * , _i , 0 , * * * , 0] P\ 


Consistent operators always have zero row-sums (J ~] ■ dj,j = 0). Thus, although of different dimension, the n — 1 
linearly independent rows of T>\ span any ?i-space stencil in the matrix T>. Represented in matrix nomenclature, this 
observation is equivalent to 

n(D\ n = niM}^ „-l[2?i]„ (181) 

where the pre- and post- subscripts are used to denote the row and column dimensions of the respective matrices. 
Substituting equation (181) into the derivative norm ||7?U||p yields 

\\W\\ 2 V = XJ T V T VVU = U t V^M t VMVi\J (182) 

We know that ||2?U||^, is bounded. If the matrix Ai T VAi is a proper norm of rank n — 1, then ||2?iU||j follows 
immediately from equation (182) by equivalence of norms. 
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By construction, the square matrix T> is of dimension n and of rank n— 1, and the matrix M. is at most rank n— 1. 
The matrix T> \ is rank n — 1, which implies from equation (182) that the matrix M. is rank n — 1. Given the full rank 
of V , the resultant matrix M T VM is symmetric positive definite, is of rank n — 1, and constitutes a proper norm. 
Thus, ||1?U||), immediately implies boundedness of ||2?iU||^-, or more precisely 

liming = c mpm (h) \\ViV\\ 2 x □ (183) 

Substituting (183) into Theorem (5) and choosing e = c mprn produces the desired pointwise estimate of the solution. 

mil < (Cmpm + r^nuf + C TOJ , m -1 ||r>U|| 2 (184) 

E Comparison of grids and geometries for polynomial collocation meth- 
ods. 

Many accuracy studies comparing liigh-order methods found in the FEM and HOFDM literature, are performed on 
uniform grids with periodic boundary conditions. While valid, the results from these studies may not be sufficiently 
general to accurately represent their behavior on nonuniform and/or nonperiodic grids. For example, the spectral 
collocation methods used in this study, often produced anomalous superconvergent behavior for uniform and/or periodic 
grids. Furthermore, the superconvergence is sometimes confined to some variables, while others converge at design 
order or less. The following studies document this complex convergence behavior for the spectral collocation methods 
used in this study. 

Consider tables 5 and 6 comparing the convergence rates of polynomial orders p = 2 to p = 5. The interface 
parameters used in the study are a = 0, /? = 1 and 7 = 0. Polynomial orders p = 2 and p = 5 superconverge at a 
rate of p + 2 for the uniform grid case, but at a rate of p + 1 for polynomial orders p = 3 and p = 4. The nonuniform 
grid case converges at the suboptimal rate of p for polynomial orders p = 2 and p = 4, at a rate of p + 1 for p = 3 and 
at p + 2 for p = 5. (In some circumstances, this anomalous behavior can be attributed to fortuitous cancellation of 
the leading order truncation terms arising in the interface derivatives used to build the penalty. Cancellation occurs 
when the truncation terms are equal. Thus, it depends on the polynomial order and does not occur on nonuniform 
grids.) Next consider a = 0, and /3 = 0 or 0 = 1. These two cases compare the nonuniform grid convergence 
behavior on the periodic and nonperiodic domains. Several of the periodic cases superconverge, while the nonperiodic 
case experiences suboptimal convergence p for all polynomial orders. 

Tables 7, and 8 show the convergence results for the four methods, four domains, and the four polynomial orders. 
Once again note the anomalous convergence behavior in table 8, exhibited for the uniform grid, nonperiodic case. This 
time for polynomials of odd order, the V integral error is superconvergent for values of f3 = j. Stranger still, is that 
this superconvergence does not improve the convergence behavior of the solution, which still converges at a suboptimal 
rate p. 

These studies makes it clear that the most demanding (and realistic) tests are those performed using nonperiodic 
boundary conditions on the nonuniform grid. Thus, in subsequent comparisons of spectral collocations methods, BC’s 
that are nonperiodic, with nonuniform grids will be the primary focus of attention. 
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Penalty 

Inflow-Outflow : Nonuniform 

Parameters 

7 = 0.000 

7 = 0.010 

7 = 0.025 

and Norms 

P = 2 

P = 3 

P = 4 

P = 5 

P = 2 

p = 3 

P = 4 

p= 5 

P = 2 

p = 3 

P = 4 

P = 5 

s 

II 

0 

II 

0 




II U - u\\ 2 

- 2.00 

-3.01 

-4.00 

-5.02 

-2.06 

-2.96 

-3.95 

-6.87 

-2.08 

-3.02 

-4.99 

-7.34 

\\U - m||oo 

- 2.01 

-2.98 

-4.00 

-5.09 

- 1.88 

-3.09 

-4.28 

-7.33 

-1.91 

-3.20 

-4.64 

-7.54 

\\U - u\\ v 

-1.93 

-3.12 

-3.97 

-5.15 

-3.35 

-2.32 

-5.10 

-5.98 

-3.36 

- 2.66 

-4.74 

-6.04 

IEHEBI 




II U - n\\ 2 

-1.97 

-3.02 

-3.97 

-5.00 

- 2.20 

-3.12 

-4.11 

-4.99 

-2.25 

-3.11 

-4.20 

-5.41 

I \U - u||oo 

- 1.88 

-3.14 

-3.99 

-5.04 

-2.09 

-3.11 

-4.00 

-5.40 

-2.30 

-3.06 

-4.37 

-5.30 

\\U - u\\-p 

-1.84 

-3.21 

-3.93 

-5.32 

-1.98 

-3.23 

-4.05 

-5.18 

-2.31 

-3.21 

-4.15 

-5.20 

a = 1 ; /? = 0 




II U - u || 2 

-3.01 

-4.00 

-4.93 

-6.03 

-2.57 

-3.58 

-4.68 

-5.55 

-2.57 

-3.60 

-4.51 

-5.49 

\\U - u||oo 

-3.01 

-3.95 

-4.92 

- 6.01 

-2.35 

-3.91 

-4.70 

-5.45 

-2.35 

-4.06 

-4.63 

-5.22 

\\U - u\\-p 

-4.04 

-4.01 

- 6.01 

-5.95 

-3.04 

-4.21 

-6.84 

-6.36 

-3.04 

-4.22 

-5.38 

-6.43 

a = 1 ; /? = j 




II U - u || 2 

-1.96 

-3.02 

-3.98 

-4.99 

-2.19 

-3.12 

-4.10 

-4.99 

-2.25 

-3.10 

-4.20 

-5.43 

\\U - u||oo 

-1.90 

-3.03 

-3.99 

-5.04 

-2.08 

-3.11 

-3.99 

-5.39 

-2.18 

-3.05 

-4.31 

-5.28 

\\U - u\\ v 

-1.84 

-3.20 

-3.93 

-5.32 

-1.97 

-3.23 

-4.04 

-5.18 

- 2.20 

-3.21 

-4.11 

-5.22 


Table 2: Convergence Study of the Linear Wave Equation comparing the effects of 7 (7 = 0). 
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Penalty 

Inflow-Outflow, Nonuniform Domain 

Parameters 

e = 1.00 

e = 0.10 

e = 0.01 

and Norms 

p = 2 

p = 3 

p = 4 

p = 5 

p = 2 

p = 3 

P — 4 

p = 5 

p = 2 

p = 3 

p = 4 

P — 5 

s 

II 

o 

II 

o 




II U - u || 2 

-1.98 

-4.01 

-4.00 

-5.84 

-1.96 

-3.98 

-4.04 

-6.01 

-2.39 

-3.74 

-4.95 

-5.92 

\\U u \\ oo 

-1.97 

-4.01 

-4.01 

-5.85 

-2.01 

-3.98 

-4.09 

-6.01 

-2.90 

-3.69 

-5.20 

-5.88 

|| vu - u x || 2 

-1.99 

-3.00 

-4.00 

-4.99 

-1.99 

-3.00 

-3.99 

-4.99 

-1.94 

-3.10 

-4.14 

-5.15 

I U - u\\-p 

-1.99 

-4.03 

-4.01 

-5.91 

-1.94 

-3.99 

-3.97 

-6.00 

-1.70 

-3.70 

-4.13 

-5.80 

a = 0;/3 = j 




II U - u || 2 

-2.01 

-3.02 

-4.02 

-5.12 

-2.98 

-3.36 

-4.85 

-5.56 

-3.38 

-4.38 

-5.20 

-6.40 

|| U - u|joo 

-2.01 

-3.02 

-4.02 

-5.08 

-2.75 

-4.09 

-4.99 

-5.96 

-2.95 

-4.47 

-4.81 

-6.54 

I T>U - u x || 2 

-2.02 

-3.00 

-4.01 

-4.99 

-2.04 

-3.02 

-4.00 

-5.00 

-2.44 

-3.33 

-4.22 

-5.26 

\\U - u\ \v 

-2.03 

-5.01 

-4.04 

-6.42 

-2.03 

-4.36 

-3.78 

-6.17 

-1.63 

-4.27 

-4.49 

-6.37 

a = l;/3 = 0 




II u - «|| 2 

-1.96 

-4.00 

-4.00 

-5.83 

-1.80 

-3.90 

-4.00 

-5.96 

-2.30 

-3.47 

-4.70 

-5.69 

\\U ^|| 00 

-1.96 

-4.00 

-4.00 

-5.84 

-1.85 

-3.90 

-4.06 

-5.95 

-2.74 

-3.41 

-4.87 

-5.64 

II vu - u x || 2 

h- 1 

io 

00 

-3.00 

-4.00 

-4.99 

-1.94 

-2.98 

-3.99 

-4.98 

-1.93 

-3.07 

-4.12 

-5.15 

\\U - u||-p 

-1.97 

-4.02 

-4.01 

-5.91 

-1.78 

-3.92 

-3.94 

-5.95 

-1.69 

-3.61 

-3.88 

-5.62 

a = l; /? = i 




lit/ - u|| 2 

-2.00 

-3.00 

-4.01 

-5.09 

-2.97 

-3.36 

-4.85 

-5.55 

-3.37 

-4.38 

-5.20 

-6.40 

- « |oo 

-2.00 

-3.01 

-4.00 

-5.04 

-2.73 

-4.08 

-4.99 

-5.96 

-2.95 

-4.47 

-4.81 

-6.54 

ijw - 2 

-2.01 

-3.00 

-4.00 

-4.99 

-2.04 

-3.02 

-4.00 

-5.00 

-2.44 

-3.33 

-4.22 

-5.26 

\\U - u\\ v 

-2.02 

-5.01 

-4.00 

-6.41 

-2.01 

-4.35 

-3.78 

-6.17 

-1.59 

-4.27 

-4.49 

-6.37 


Table 3: Convergence Study: Effects of e in the linear Burgers’ equation, using the primal (B-O) formulation. 
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Penalty 

Inflow-Outflow, Nonuniform Domain 

Parameters 

e = 1.00 

e = 0.10 

e = 0.01 

and Norms 

p = 2 

p = 3 

p = 4 

p = 5 

p = 2 

p = 3 

p = 4 

p = 5 

p = 2 

p = 3 

p = 4 

P — 5 

s 

II 

o 

II 

o 




II U - u|| 2 

-2.97 

-4.07 

-4.96 

-5.84 

-2.88 

-3.89 

-4.94 

-5.92 

-3.42 

-4.24 

-5.37 

-6.10 

\\U u \\ oo 

-2.96 

-3.99 

-4.99 

-5.85 

-2.59 

-3.85 

-4.75 

-5.87 

-3.43 

-4.12 

-5.22 

-5.86 

|| vu - u x || 2 

-2.01 

-3.03 

-4.02 

-5.08 

-2.02 

-3.03 

-4.03 

-5.08 

-2.28 

-3.41 

-4.35 

-5.45 

I U - u\\-p 

-2.81 

-3.81 

-5.86 

-7.28 

00 

T— 1 

-3.94 

-3.72 

-5.91 

-2.44 

-3.91 

-3.19 

-6.00 

a = 0;/3 = j 




II U - u|| 2 

-2.90 

-3.96 

-4.91 

-5.82 

-3.01 

-4.04 

-4.98 

-6.00 

-3.72 

-4.39 

-5.56 

-6.30 

|| U - u|joo 

-2.91 

-3.92 

-5.02 

-5.90 

-3.43 

-4.04 

-4.99 

-6.02 

-3.41 

-4.42 

-5.54 

-6.28 

I T>U - u x || 2 

-2.00 

-3.00 

-3.98 

-5.00 

-1.99 

-2.99 

-3.98 

-5.02 

-2.59 

-3.39 

-4.44 

-5.45 

\\U - u\ \v 

-2.61 

-3.84 

-5.62 

-6.50 

-0.33 

-4.21 

-4.57 

-6.08 

-2.42 

-4.25 

-4.10 

-6.51 

a = l;/3 = 0 




||C/ - «|| 2 

-2.95 

-4.06 

-4.95 

-5.84 

-2.81 

-3.82 

-4.91 

-5.88 

-3.18 

-3.99 

-5.10 

-5.67 

\\U ^|| 00 

-2.93 

-3.98 

-4.98 

-5.84 

-2.46 

-3.78 

-4.69 

-5.83 

-3.33 

-3.89 

-4.98 

-5.67 

II vu - u x || 2 

-2.01 

-3.03 

-4.02 

-5.08 

-2.01 

-3.02 

-4.01 

-5.05 

-2.01 

-3.12 

-4.09 

-5.13 

\\U - u||-p 

-2.79 

-3.79 

-5.86 

-7.30 

-1.98 

-3.88 

-3.72 

-5.87 

-2.78 

-3.70 

-3.51 

-5.78 

a = l; /? = i 




||t/ - u|| 2 

-2.90 

-3.94 

-4.90 

-5.80 

-3.00 

-4.04 

-4.98 

-6.00 

-3.72 

-4.39 

-5.56 

-6.30 

II t/ - w joo 

-2.92 

-3.90 

-5.02 

-5.88 

-3.41 

-4.03 

-4.98 

-6.02 

-3.41 

-4.42 

-5.54 

-6.28 

|T>t/ - u x || 2 

-1.99 

-2.99 

-3.98 

-5.00 

-1.99 

-2.99 

-3.98 

-5.02 

-2.59 

-3.39 

-4.44 

-5.45 

e/ - u||-p 

-2.59 

-3.79 

-5.63 

-6.52 

2.18 

-4.20 

-4.56 

-6.07 

-2.41 

-4.25 

-4.10 

-6.50 


Table 4: Convergence Study: Effects of e in the linear Burgers’ equation, using the flux (LDG) formulation. 
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Penalty 


Parameters 


and Norms 


Periodic Domain 


Uniform Grid 


p=2p=3p=4p=5 


Nonuniform Grid 


p=2p=3p=4p=5 


a = 

0 

II 

0 

\\U 

- u h 

\\U 

^||oo 

\\U 

- u\\-p 

a = 

WBM 

II U 

- U || 2 

\\U 

^11 00 

\\U 

- u\\ v 

a = 

i;/? = o 

\\U 

- u\\ 2 

\\U 

^11 00 

\\U 

— u\\j> 


-4.05 

-4.00 

-4.96 

-7.07 

-3.99 

-4.03 

-5.30 

-6.67 

-3.99 

-4.00 

-5.99 

- 6.02 


-1.97 

-4.30 

-3.93 

-7.00 

-1.95 

-4.20 

-3.96 

-6.70 

-4.00 

-4.00 

-5.99 

-5.98 


-3.00 

-4.01 

-5.07 

- 6.00 

-3.00 

-4.01 

-5.06 

-5.99 

-4.03 

-4.00 

-5.99 

-5.95 


-3.01 

-4.00 

-4.98 

-6.04 

-3.01 

-3.98 

-4.94 

-6.03 

-4.03 

-4.00 

- 6.00 

-5.97 


II U - u || 2 

-4.00 

-4.25 

-4.90 

-7.05 

-1.96 

-4.60 

-3.98 

-6.07 

| \u - It |oo 

-4.00 

-4.06 

-5.12 

-6.67 

-1.97 

-4.23 

-3.96 

-6.65 

\\U - ujj-p 

-4.00 

-4.00 

-5.99 

- 6.01 

-4.00 

-4.00 

- 6.00 

-5.98 


Table 5: Convergence Study: Scalar Wave Equation on a uniform grid (7 = 0). 


Penalty 

Inflow- Outflow Domain 

Parameters 

Uniform Grid 

Nonuniform Grid 

and Norms 

P= 2 

p = 3 

p = 4 

P = 5 

p = 2 

p = 3 

p = 4 

P= 5 

s 

II 

0 

II 

0 



II U - w || 2 

-2.95 

-2.99 

-4.99 

-5.01 

- 2.00 

-3.01 

-4.00 

-5.02 

| \U - u||oo 

-3.03 

-3.01 

-4.95 

-5.04 

- 2.01 

-2.98 

-4.00 

-5.09 

\\U - u\\ v 

-3.13 

-3.11 

-4.95 

-5.02 

-1.93 

-3.12 

-3.97 

-5.15 




II u - «|| 2 

-3.04 

-3.03 

-4.98 

-4.98 

-1.97 

-3.02 

-3.97 

-5.00 

II U ^|| 00 

-3.15 

-3.03 

-5.01 

-4.84 

- 1.88 

-3.14 

-3.99 

-5.04 

II U - u\\ v 

-3.09 

-3.16 

-4.97 

-5.11 

-1.84 

-3.21 

-3.93 

-5.32 

a = l;/3 = 0 



II U - u || 2 

-2.99 

-4.01 

-4.98 

- 6.00 

-3.01 

-4.00 

-4.93 

-6.03 

| \U - ulloo 

-2.99 

-4.01 

-4.99 

-5.98 

-3.01 

-3.95 

-4.92 

- 6.01 

\\U - u\\v 

-4.02 

-4.00 

-5.98 

-5.92 

-4.04 

-4.01 

- 6.01 

-5.95 

a = l;/3= j 



II u - «|| 2 

-3.02 

-3.00 

-4.98 

-4.97 

-1.96 

-3.02 

-3.98 

-4.99 

II U ^|| 00 

-3.01 

-3.03 

-5.00 

-4.88 

-1.90 

-3.03 

-3.99 

-5.04 

\\U - u\\ v 

-3.09 

-3.11 

-4.97 

-5.17 

-1.84 

-3.20 

-3.93 

-5.32 


Table 6 : Convergence Study: Scalar Wave Equation on a nonuniform grid (7 = 0). 
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Penalty 

Periodic Domain 

Parameters 

Uniform Grid 

Nonuniform Grid 

and Norms 

p = 2 

p = 3 

p = 4 

P = 5 

p = 2 

P = 3 

p = 4 

p = 5 

S 

II 

0 

II 

0 



II u - «|| 2 

- 1.99 

- 4.00 

- 4.02 

- 6.40 

- 1.99 

- 4.00 

- 4.01 

- 5.82 

|| u - u||oo 

- 1.98 

- 4.00 

- 4.02 

- 6.37 

- 1.99 

- 4.00 

- 4.01 

- 5.90 

II VU - u x || 2 

- 1.99 

- 3.00 

- 4.01 

- 5.00 

- 1.99 

- 3.00 

- 4.00 

- 5.00 

\\U - u\\ v 

- 1.98 

- 4.00 

- 4.02 

- 6.43 

- 1.98 

- 4.00 

- 4.01 

- 5.88 

a = ();0=j 



II U - u\\ 2 

- 2.00 

- 3.01 

- 4.03 

- 5.05 

- 2.00 

- 3.00 

- 4.01 

- 4.99 

| \u - u||oo 

- 2.00 

- 3.01 

- 4.03 

- 5.06 

- 2.00 

- 3.01 

- 4.01 

- 5.05 

|| VU - u x || 2 

- 2.01 

- 3.00 

- 4.01 

- 5.00 

- 2.01 

- 3.00 

- 4.00 

- 5.00 

\\U - u\ \ v 

- 2.00 

- 4.00 

- 4.03 

- 6.54 

- 2.00 

- 4.00 

- 4.01 

- 5.88 

a = 1 ; (3 = 0 



II U - u\\ 2 

- 1.97 

- 4.00 

- 4.01 

- 6.40 

- 1.97 

- 4.00 

- 4.00 

- 5.81 

|| u - u||oo 

- 1.97 

- 4.00 

- 4.01 

- 6.37 

- 1.97 

- 4.00 

- 4.00 

- 5.89 

|| VU - U x || 2 

- 1.97 

- 3.00 

- 4.01 

- 5.00 

- 1.98 

- 3.00 

- 4.00 

- 4.99 

\\U - u\\ v 

- 1.96 

- 4.00 

- 4.01 

- 6.43 

- 1.97 

- 4.00 

- 4.00 

- 5.88 

a = i; /3 = i 



II U - u\\ 2 

- 1.99 

- 3.00 

- 4.03 

- 5.05 

- 1.99 

- 3.00 

- 4.01 

- 4.99 

| \u - u|joo 

- 1.99 

- 3.01 

- 4.03 

- 5.06 

- 1.98 

- 3.01 

- 4.01 

- 5.04 

II VU - u x || 2 

- 1.99 

- 3.00 

- 4.01 

- 5.00 

- 2.00 

- 3.00 

- 4.00 

- 4.99 

\\U - u\\ v 

- 1.99 

- 4.00 

- 4.03 

- 6.54 

- 1.99 

- 4.00 

- 4.01 

- 5.88 


Table 7: Convergence Study of the Linear Burgers’ equation using the Baumann-Oden scheme on a periodic domain 
(7 = 0 ). 
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Penalty 

NonPeriodic Domain 

Parameters 

Uniform Grid 

Nonuniform Grid 

and Norms 

P = 2 

p = 3 

p = 4 

P = 5 

p = 2 

P = 3 

p = 4 

P = 5 

S 

II 

0 

II 

0 



II u - u || 2 

-1.97 

-4.01 

-4.01 

-6.18 

h - 1 

io 

00 

-4.01 

-4.00 

-5.84 

|| u - u||oo 

-1.97 

-4.02 

-4.01 

-6.20 

-1.97 

-4.01 

-4.01 

-5.85 

II vu - Ux || 2 

-1.99 

-3.00 

-4.01 

-5.00 

-1.99 

-3.00 

-4.00 

-4.99 

\\U - u\ \v 

-1.98 

-4.02 

-4.02 

-6.40 

-1.99 

-4.03 

-4.01 

-5.91 

wasmsm 



II U - u || 2 

-2.01 

-2.95 

-4.02 

-4.99 

-2.01 

-3.02 

-4.02 

-5.12 

I \U - m||oo 

-2.00 

-3.02 

-4.03 

-5.07 

-2.01 

-3.02 

-4.02 

-5.08 

|| VU - Uxh 

-2.02 

-3.00 

-4.01 

-5.00 

-2.02 

-3.00 

-4.01 

-4.99 

\\U - u\ \ v 

-2.03 

-6.02 

-4.02 

-7.45 

-2.03 

-5.01 

-4.04 

-6.42 

a = 1; (3 = 0 



II u - «|| 2 

-1.96 

-4.01 

-4.01 

-6.18 

-1.96 

-4.00 

-4.00 

-5.83 

|| u - w||oo 

-1.95 

-4.01 

-4.01 

-6.19 

-1.96 

-4.00 

-4.00 

-5.84 

II vu - « x || 2 

-1.97 

-3.00 

-4.01 

-5.00 

-1.98 

-3.00 

-4.00 

-4.99 

\\U - u||-p 

-1.96 

-4.01 

-4.01 

-6.40 

-1.97 

-4.02 

-4.01 

-5.91 

a = 1; /3 = J 



II U - «|| 2 

-1.99 

-2.94 

-4.02 

-4.99 

-2.00 

-3.00 

-4.01 

-5.09 

I \u - u|joo 

-1.99 

-3.02 

-4.03 

-5.07 

-2.00 

-3.01 

-4.00 

-5.04 

II vu - « x || 2 

-2.01 

-3.00 

-4.01 

-5.00 

-2.01 

-3.00 

-4.00 

-4.99 

\\U - u\\ v 

-2.01 

-6.07 

-4.04 

-7.46 

-2.02 

-5.01 

-4.00 

-6.41 


Table 8: Convergence Study of the Linear Burgers’ equation using the Baumann-Oden scheme on a nonperiodic domain 
(7 = 0 ). 


47 





































































































































































Wave Equation 



Figure 5: Comparison of 4 th - and 6 ^ -order finite difference schemes on wave equation. 
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ELEMENTS 


Figure 7: Comparison of Baumann-Oden and LDG schemes for multiple penalty parameters on the linear Burgers’ 
equation: P = 2. 



Figure 8: Comparison of Baumann-Oden and LDG schemes for multiple penalty parameters on the linear Burgers’ 
equation: P = 3. 
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30 60 90 120 

ELEMENTS 


Figure 9: Comparison of Baumann-Oden and LDG schemes for multiple penalty parameters on the linear Burgers’ 
equation: P = 4. 



ELEMENTS 


Figure 10: Comparison of Baumann-Oden and LDG schemes for multiple penalty parameters on the linear Burgers’ 
equation: P = 5. 
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Figure 11: Comparison of Baumann-Oden and LDG schemes for multiple polynomial orders on the linear Burgers’ 
equation: a = 1, /? = 0. 
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